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Natural and artificial light harvesting processes have recently gained new interest. Signatures of long lasting coherence 
in spectroscopic signals of biological systems have been repeatedly observed, albeit their origin is a matter of ongoing 
debate, as it is unclear how the loss of coherence due to interaction with the noisy environments in such systems is 
averted. Here we report experimental and theoretical verihcation of coherent exciton-vibrational (vibronic) coupling 
as the origin of long-lasting coherence in an artihcial light harvester, a molecular J-aggregate. In this macroscopically 
aligned tubular system, polarization controlled 2D spectroscopy delivers an uncongested and specihc optical response 
as an ideal foundation for an in-depth theoretical description. We derive analytical expressions that show under which 
general conditions vibronic coupling leads to prolonged excited-state coherence. 


Introduction 

The remarkably high efficiency in photosynthesis, where 
nine out of ten absorbed photons reach the reaction center, is 
a fascinating held of modern research. In such photosynthetic 
complexes, structure, dynamics and function are inextricably 
linked. A conserved building block comprises strongly ab¬ 
sorbing pigments arranged in close proximity to one another 
by a surrounding protein scaffoldi*^. Typical inter-pigment 
distances are of order of 10 A and photon absorption leads 
to the formation of delocalized excited electronic states (exci- 
tons) shared by two or more pigment molecules. Exciton cre¬ 
ation, migration and trapping are central to the functionality 
of a photosynthetic apparatus. The controlled and adjustable 
arrangement of the pigments tunes the electronic network and 
the properties of its interaction with the vibrational environ¬ 
ment that is associated with either the pigments or the protein. 
The detailed balance of these properties determines the effi¬ 
ciency of light harvesting systems^. 

Exciton dynamics can be efficiently probed by two- 
dimensional (2D) electronic spectroscopy^. This technique 
revealed oscillatory signals in the spectral response of a wide 
variety of photosynthetic aggregates^. Initially ascribed to 
excitonic beatings, oscillations have been found to persist up 
to several hundreds of femtoseconds at room temperature^J^. 
This time scale exceeds typical dephasing rates in the con¬ 
densed phase and becomes comparable to exciton transfer 
timesi, thus posing the question of the nature and functional 
relevance of these coherences^. Unfortunately, the complex 
structure of 2D signals makes the unambiguous identihcation 
of the underlying mechanisms that support such long-lived co¬ 
herences a challenging task and several hypotheses to explain 
them have been formulatedii.^. The different approaches can 
be classihed into theories including coherent interaction of ex- 
citons with intra-pigment vibrationsiivl.^ and theories focus¬ 
ing on incoherent exciton-protein interaction such as corre¬ 
lated fluctuationsi^J^. It is possible that some of these mech¬ 
anisms may coexist on certain time scales and that one or an¬ 


other may become dominant depending on the system under 
consideration. 

In this work, we show that the relatively simple excitonic 
structure of a molecular J-aggregate provides an ideal test 
case to identify the microscopic mechanism behind long- 
lived oscillations in electronic 2D-signals. The investigated J- 
aggregate is tubular and aligns along the sample’s flow direc¬ 
tion when in solution. Additionally, the J-aggregate exhibits 
excitonic bands with roughly orthogonal transition dipole mo¬ 
ments. It is this combination of perpendicular excitonic tran¬ 
sitions and macroscopic alignment that makes electronic 2D- 
spectroscopy with polarization-controlled excitation pulses an 
ideal tool to study coherence effects between the excitonic 
bands. This approach signihcantly reduces the complexity of 
retrieved 2D signals, leading to only two peaks with oscilla¬ 
tory components in specihc regions of the 2D-maps, i.e. one 
on the diagonal and one as a cross-peak for non-rephasing and 
rephasing signal components, respectively. Employing a vi¬ 
bronic model, we derive analytical expressions that show how 
system parameters such as electronic decoherence rates and 
exciton-vibrational resonance determine the amplitude and 
lifetime of oscillatory signals. Eitting the analytical expres¬ 
sions to measured data, the vibronic model achieves quanti¬ 
tative agreement with experimental observations. Concerning 
potential functional relevance of the observed oscillations, we 
show that the long-lived oscillatory signals in our system are 
dominated by excited-state coherence rather than ground-state 
coherence. 

Results 

The system. J-aggregates of cyanine dyes are promis¬ 
ing candidates for artihcial antenna systems22;i2^. They are 
chemically versatile and self-assemble into various extended 
supramolecular structures in aqueous solution^!. Here a sys¬ 
tem that can be considered a macroscopically aligned syn¬ 
thetic light harvester was studied, namely a molecular J- 
aggregate of C803-monomers whose aggregation behavior 
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FIG. 1. C803 and polarization controlled 2D spectroscopy, a. 

Wire-guided window-free jet used for sample circulation, along with 
a schematic of the double-layered structure of the C803-aggregate. 
The aggregates align along the flow direction (white arrow). The 
transition dipole directions of bands 1-3 are displayed by arrows, 
which are mainly polarized along the tube axis (bands 1 and 2 shown 
in blue) or perpendicular to the axis (band 3 shown in orange), b. Ab¬ 
sorption spectra with light polarized parallel (blue) and perpendicular 
(orange) to the flow direction, c. Non-resonant Raman spectra of the 
C803-monomer (black line) and aggregate (grey area). The vibra¬ 
tional frequencies Vi and V 2 are close to the exciton energy splitting 
between bands 1 and 3 and bands 2 and 3, respectively, d. Polariza¬ 
tion controlled 2D spectroscopy with three excitation pulses (ki to 
k^) and a local oscillator (LO) for heterodyne detection of the signal 
field, depicted as an oscillating line. Polarization orientation (0° or 
90°) is given with respect to the longitudinal axis of aligned C803. 


is well known22i22. As revealed by cryogenic transmission 
electron microscopy^, the aggregate structure is best de¬ 
scribed as a double-layered nanotube with outer diameter 
~ 11 nm and lamellar spacing of ~ 2.2 nm between the chro- 
mophore layers. Additionally, superhelical bundles of these 
tubes can also form, though the addition of polyvinyl alcohol 
(PVA) inhibits this process and thereby avoids single-layered 
tube formation^! and maintains a stable solution over several 
weeks^. A drawing of the J-aggregate under investigation, 
from here on referred to as C803, is shown in Fig. la. The 
bilayer configuration of C803 allows the effect of different 
decoherence rates to be studied as the outer solvent-exposed 
layer shows faster decoherence than the inner protected layer. 

The structural properties of the aggregate are remarkable; 
the 11 nm outer diameter is contrasted by a length of several 
micrometers. Circulating solvated C803 with a wire-guided 
jet (Fig.la) leads to a macroscopic orientation of the tubes 


because the longitudinal axis preferentially aligns along the 
flow direction. This creates anisotropy for linearly polarized 
light, as shown in Fig. lb. Linear dichroism measurements^ 
and redox-chemistry studies^ assign bands 1 and 2 to longi¬ 
tudinal transitions localized upon the inner and outer cylin¬ 
ders, respectively (Fig. la). Transitions to band 3 are pref¬ 
erentially polarized perpendicular to the long axis of C803 
and are shared by both layers. A detailed description of sam¬ 
ple preparation methods and band assignments is given in the 
Supplementary Notes 1 and 2. 

Fitting the well-defined absorption peaks of C803 with 
Lorentzian functions (see Supplementary Note 2) reveals an 
exciton energy difference between bands 1 and 3 of AQ 31 « 
690 cm“^ and AQ 32 ~ 460 cm * for bands 2 and 3. Both 
exciton energy splittings are close to vibrational frequencies 
vi » 668 cm“^ and V 2 ~ 470 cm“* observed in non-resonant 
Raman spectra^ (Fig.lc). These vibrational frequencies are 
measured in both the monomer and aggregate Raman spectra, 
i.e. they are not aggregation induced Raman bands. Strongly 
enhanced modes at similar energies were observed in resonant 
Raman spectra of a related cyanine dye, and can be assigned to 
out-of-plane vibrations^. Such out-of-plane vibrations were 
shown to couple strongly to excitons^. The quasi-resonance 
between the vibrational frequencies vi and V 2 and exciton en¬ 
ergy splittings Af 23 i and Afl 32 provides us with an interesting 
scenario of possible coherent interaction between bands (exci- 
tons) and vibration s*^'^^'*"^'^''^^ . Such exciton-vibrational cou¬ 
pling induces vibronioi^ and vibrational coherences'^, which 
can both lead to long-lived beating signals in 2D spectra. 
Here we emphasize that coherence in the electronic excited- 
state manifold is referred to as vibronic and in the ground- 
state manifold as vibrational. Identifying the dominant con¬ 
tribution is of fundamental importance because only vibronic 
coherence, which manifests in excited state dynamics, can 
enhance exciton transport and thus support light-harvesting 
function^^— . 

Experimental results. The absorption spectrum of a light 
harvesting system may be heavily congested because of over¬ 
lapping excitonic bands and the resulting 2D-signal would ex¬ 
hibit significant overlap between diagonal and cross peaks, 
thereby impeding further analysis. It has been suggested to 
employ laser pulses of different relative polarization to selec¬ 
tively address relevant excitation pathways to obtain a clearer 
2D signal^. However, the advantage of polarization con¬ 
trolled 2D spectroscopy has been limited by the isotropic na¬ 
ture of the investigated samples (an ensemble). In the experi¬ 
ment presented here, these problems are circumvented by the 
measurement of the macroscopically aligned C803. The tran¬ 
sition dipole moments of bands 1 and 2 are preferentially par¬ 
allel to the longitudinal axis while band 3 is orthogonal, thus 
allowing for optimal polarization selectivity. This combina¬ 
tion reduces the obtained 2D maps to only two relevant peaks 
with negligible overlap and an up to 30 times stronger signal 
intensity as compared to the isotropic caseii. 

The ideal pulse sequence to isolate beating signals between 
states with orthogonal transition dipole moments, i.e. bands 
1 and 3 in the present case, is depicted in Fig. Id, where 
the phase-matched direction for measuring rephasing spectra 
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FIG. 2. Experimental and theoretical 2D spectra, a, b. The 

Fourier-transform amplitude maps of non-rephasing and rephasing 
spectra at W 2 = 705 ± 20cm“*, which reveal the presence of a non¬ 
rephasing diagonal peak Nil and a rephasing cross-peak R31. These 
peaks stem from the coherent interaction of bands 1 and 3 with the 
quasi-resonant vibrational mode with frequency vi a; 668 cm“*. The 
amplitude of Nil is about three times larger than R31. The line- 
shape of Nil is symmetric along both Wi- and Ws-axes, while that 
of R31 is elongated along tui-axis. c, d, The simulated spectra at 
W 2 = 705cm * with Nil and R31. e, The FT amplitude map at 
W 2 = 462±20 cm“* reveals coherent interaction of bands 2 and 3 with 
the quasi-resonant vibrational mode with frequency V 2 ~ 470 cm“'. 
However, as depicted in f, the associated non-rephasing peak N22 
at 0)13 ~ 16670cm“' is weak and only amounts to 5% of Nil at 
0)2 = 705 ± 20 cm“* (see a). The diagonal peak at 0)13 ~ 16400 cm“* 
in e stems from Nil, with a peak centered at 0)2 = 705±20cm“*, but 
broad enough to appear at 0)2 = 462 ± 20cm“'. All measurements 
were carried out at room temperature. 


is displayed: non-rephasing spectra can be measured along 
the same phase-matched signal direction by changing the or¬ 
der of the first two pulses (see Methods). After subtraction 
of the non-oscillatory background, we performed a Fourier 
transformation along waiting time (2 for all points on the 
two-dimensional ((x)i,W 3 )-map. The resulting (x) 2 -plots al¬ 
low the lineshape of beating signal with frequency C 02 to be 
visualized as a function of position in (wi, W 3 )-space. The 
slice at the exciton energy splitting between bands 1 and 3 
(ai 2 = 705 + 20cm“' with the experimental resolution of 
+20cm ’) reveals a non-rephasing diagonal peak Nil and a 
rephasing cross-peak R31 as shown in Figs. 2a and b, respec¬ 
tively. N11 is centered at (wi, W 3 ) = (Qi, Qi) with exciton en¬ 
ergy Qi 16405 cm ' of band 1 and a symmetric linewidth 
2Fji ai 130cm“' along both oJi- and (U 3 -axes (Fig. 2a). The 
center of R31 is located at (<^ 1 , 0 ) 3 ) = (Q 3 ,f 2 i) with ex¬ 
citon energy ^ 17125 cm“* of band 3 and asymmetric 
linewidths 2 Fg 3 300cm“' and 2Tg\ » 130cm“' along wi- 
and a) 3 -axes, respectively (Fig. 2b). In peak amplitude, R31 
is approximately 30 % of N11. Turning to the W 2 -slice cor¬ 
responding to the energy splitting between bands 2 and 3, 
(aJ 2 = 462 + 20cm“'), Figs.2e and f reveal a diagonal non¬ 
rephasing peak N22, which is centered at (wi, 013 ) = (Q 2 , ^ 2 ) 
with the exciton energy Q 2 ~ 16672 cm ' of band 2 and a 
symmetric linewidth 2 Fg 2 ~ 225 cm“' along a>i- and W 3 -axes. 
The amplitude of N22 is only 5 % of N11. 

Theoretical model. In order to describe the long-lived os¬ 
cillations in N11 and R31, a vibronic model is employed that 
describes the coupling of bands 1 and 3 to a quasi-resonant 
vibrational mode with frequency vi. Consider a system with 
electronic ground state and excited states for bands 1 and 
3, denoted by |l,i) and |3j;), respectively, where k - Q and 1 
denote the vibrational ground and excited state, respectively 
(Fig. 3a). The vibronic coupling between the quasi-resonant 
states |3o) and |li) leads to unnormalized vibronic eigenstates 
(3o| = (3ol + e{lil and (ii| = {li| - e{3o|. Here, e represents 
the degree of vibronic mixing defined by 


e = /vi V^(iAvi -F 13 ) 

where Avi = (Q 3 - Qi) - vi denotes the detuning between 
|3o) and |li), i.e. between the exciton energy splitting and vi¬ 
brational frequency, and S \ denotes the Huang-Rhys factor of 
the vibrational mode, which in turn quantifies the strength of 
the vibronic coupling (see Supplementary Note 2 for details 
of the derivation). The electronic decoherence rate F^^ de¬ 
scribes the exponential decay rate of the coherence between 
electronic ground state and band k, while F 13 represents the 
overall exponential decay rate of the inter-exciton coherence 
between bands 1 and 3. In our model, we do not consider 
inhomogeneous broadening, which is justified by the obser¬ 
vation that the experimentally measured absorption spectrum 
is well matched to a sum of Lorentzian functions with the 
linewidths 2Tgi, (see Supplementary Note 2). This is valid 
when homogeneous broadening dominates the linewidths and 
the Huang-Rhys factors are sufficiently small, as is the case 
here. In addition, the lineshape of Nil (Fig. 2a) is not elon¬ 
gated along the diagonal u)\ = 023 , implying our 2D signal is 
dominated by homogeneous broadening. The same conclu- 
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sion is reached from analyzing 2D correlation spectra^. 

In nonlinear spectroscopy, the molecular response to laser 
excitation is described by response functions^. According to 
the vibronic model described above, the response function for 
the oscillatory signals in N11 reads 

with yUi and yU 3 denoting the transition dipole moment of bands 
1 and 3, respectively. The prefactor stems from the line- 
shape of N 11 , yv denotes the dissipation rate of the vibrations 
and 6 (i) stands for the frequency shift of the vibronic eigen¬ 
states (3ol and (li| relative to the uncoupled states {3ol and 
{111 due to the vibronic coupling (see Fig. 3a and Supplemen¬ 
tary Note 2 for further details). The coupling was found to be 
sufficiently strong to induce non-negligible vibronic mixing 
|ep =» 0.03, which leads to a long-lived beating signal in N11 
up to t 2 ~ 800 fs, as shown in Fig.3b. These results imply 
that the initial excitonic part of |lo)(3ol decays rapidly with 
1 /e decay time of F/j 66 fs, while the vibronic coherence 
|lo){iil explains a long-lived oscillatory signal in Nil: here 
|lo){3ol (|lo){lil) represent coherence between two vibronic 
states |lo) and (3ol (|lo) and {li|), respectively. 

The response function for the oscillatory contributions to 
R31 is given by 

where r^ 3 *r^i* derives from the asymmetric lineshape of R31 
(see Figs.2b and d). Here rje and rjg represent the contribution 
of excited-state vibronic coherence |lo){lil and ground-state 
vibrational coherence |go) {.gil, respectively, to the long-lived 
beating signal in R31 (see Supplementary Note 2). The vibra¬ 
tional coherence in the electronic ground-state manifold does 
not play a role in exciton transfer dynamics, but nonetheless 
modulates the 2D spectra. A fit of model parameters to ex¬ 
perimental results (Fig. 3c) shows that Irid « 2.5 | 77 ^|. This 
means the long-lived beating signal in R31 is dominated by 
the excited-state coherence |lo)(lil- The short-lived beating 
signal in R31 is induced by |lo) {3o|, as is the case forNl 1. We 
note that the signal at N 11 , with approximately three times the 
amplitude of R31, is exclusively determined by excited-state 
contributions. Details of this vibronic model and the corre¬ 
sponding Feynman diagrams for the spectral components Nil 
and R31 are discussed in the Supplementary Note 2. 

These results demonstrate how an excitonic system within 
a noisy environment can exhibit long-lasting coherent fea¬ 
tures: the observed long-lived oscillations are the result of 
coherent interaction of excitonic bands with an underdamped, 
quasi-resonant vibration. This vibronic mechanism requires 
the vibrational dissipation rate to be much slower than the 
electronic decoherence rate F 13 , which is the case for C803, 
where < (lps)“* and F 13 « ( 66 fs) *. The difference in 
electronic and vibrational decoherence rates can be rational¬ 
ized from the fact that excitons and vibrations are related to 
the motion of electrons and nuclei, respectively. The lower 
mass of electrons as compared to nuclei makes excitons more 
mobile and therefore more sensitive to environmental fluctu¬ 
ations, such as local electric fields, than vibrations. We note 


that the vibronic mixing leading to long-lived beating signals 
in 2D-ES is described by a vibronic coupling that induces co¬ 
herent energy exchange between excitons and quasi-resonant 
vibrations (see Supplementary Note 2 for further details): 

= V^(|3o)<li| + |li><3ol). 

This implies that the vibronic coupling not only induces long- 
lasting electronic excited-state coherences, but also can medi¬ 
ate population transfer between excitonic bands. In a combi¬ 
nation with thermal relaxation of exciton populations, the vi¬ 
bronic coupling may further enhance exciton population trans¬ 
fer and as a result could, in principle, have functional rele¬ 
vance in exciton transpor t*'*’^^^"*^^ . 

Interestingly, the different decoherence rates Fg 3 2Fgi of 
bands 1 and 3 lead to different amplitudes of the short-lived 
beating signals in N11 and R31 (Figs.3b and c), which are de¬ 
termined by the prefactors r~^ and F^jF^*, respectively. The 
lower decoherence rate of band 1 can be explained by band 1 
being localized on the inner layer, while band 3 is delocalized 
over both the inner and outer layersi^. As shown by the re¬ 
sponse functions for Nil and R31, the overall strength of the 
beating signals is proportional to the inverse of the electronic 
decoherence rates. It is therefore expected that the beating 
signal amplitude would diminish with an increase of the de¬ 
coherence rate. This is the case for N22, where the physical 
situation in terms of exciton-vibrational resonance (Afl 32 ~ 
V 2 ~ 470 cm“*) is equivalent to Nil (AD 31 ^ vi s» 668 cm“*). 
The crucial difference is that band 2 has a higher decoher¬ 
ence rate than band 1 , as band 2 is localized on the outer layer 
exposed to solvent^. This explains the broader linewidth of 
band 2 in absorption and 2D spectra. Using an estimated value 
ofFg 2 ~ (47 fs)“', the presented theory predicts the strength of 
N22 to be 5 % of N11 (see Supplementary Note 2), which is in 
line with the experimental observations (Fig.2f). These results 
indicate that the experimentally observed long-lived beating 
signals, induced by vibronic mixing, require adequately low 
electronic decoherence rates, highlighting that resonance be¬ 
tween exciton energy splitting and vibrational frequency alone 
is not sufficient^. 

The presented vibronic model achieves quantitative agree¬ 
ment with the experimental observations. Crucially, the con¬ 
straints imposed by the observed asymmetric decoherence 
rates rg 3 a; 2 Fgi and fast relaxation of exciton population 
in C803 on sub-picosecond timescales^ rule out incoherent 
models, where long-lived oscillations are sustained by Marko¬ 
vian correlated fluctuations (see Supplementary Note 3 for a 
detailed analysis). This further supports our conclusion that 
the observed experimental data provide evidence for vibronic 
mixing being the mechanism at play in our system. 

We note that our results do not imply that correlated fluctua¬ 
tions can be universally ruled out, as this mechanism could be 
in place in certain pigment-protein complexes. The notion of 
correlated fluctuations has been developed for photosynthetic 
complexes where pigments are embedded in a protein scaf¬ 
fold. The protein has been considered as the potential source 
of correlated fluctuations in natural light harvester s*^'*^ . For 
C803, a structural frame such as a protein scaffold is absent 
and therefore correlated fluctuations are unlikely to induce 
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FIG. 3. Vibronic model, a, We consider a vibronic model for bands 
1 and 3 coupled to a vibrational mode with frequency vi ~ 668 cm“* 
(see Supplementary Note 2). The vihronic states |A:o) and |^:i) denote 
the vibrational ground and first excited state of an electronic state 
1^:), respectively, with the single index states |g), |1) and |3) denot¬ 
ing the electronic ground state and bands 1 and 3, respectively. The 
exciton energy splitting Aflsi = fls - between bands 1 and 3 is 
quasi-resonant with the vibrational frequency V[, where the detuning 
is denoted by Avi = AQsi - vi. The exciton-vibrational coupling be¬ 
tween uncoupled states |3o) and |li) leads to vibronic eigenstates |3o) 
and 111 ), each of which is a superposition of |3o) and |li), leading to 
an energy-level shifting by Scl). b. The time trace of N11 where the 
experimental results are shown as light red circles, and the theoretical 
simulation is shown as a full red line, c. The time trace of R3 1 where 
the experimental results are shown as light blue circles, and the sim¬ 
ulated data are depicted as a full blue line. The root-mean-square 
deviation (RMSD) between the experimental results and theoretical 
simulation in b and c is 0.92 and 0.59, respectively. 


long-lived oscillatory 2D-signals, which is in line with our 
observations. 

Discussion 

We have verified, theoretically and experimentally, that co¬ 
herent vibronic coupling in the electronic excited-state mani¬ 
fold is responsible for the long-lived beating signals observed 
in 2D spectra of an artificial light harvester. The relatively 
simple electronic and vibrational structure of the investigated 
molecular aggregate along with its macroscopic alignment al¬ 
lowed us to rule out the presence of correlated fluctuations. 
The specific geometry of our system allowed us to gain fur¬ 
ther insights by illustrating the conditions under which intra¬ 
pigment vibrations can prolong electronic coherent effects. 
The moderately low decoherence rate of band 1, localized on 


the inner layer and protected from solvent, is the basis for 
exciton-vibrational coupling as the source of long-lived beat¬ 
ing signals. The outer band 2, even though resonantly coupled 
to a vibration, exhibits a higher decoherence rate and therefore 
fails to produce observable oscillations. We conclude that the 
mere resonance between excitons and vibrations does not suf¬ 
fice to explain long-lived beating signals. An adequately low 
electronic decoherence rate, determined by the interaction be¬ 
tween system and bath, is an equally important prerequisite. 

The influence of vibronic coupling on energy transport 
in molecular aggregates has been extensively studied in the 
past, as recently reviewed^. The vibronic coupling has re¬ 
cently gained new interest (see ref4^ for a recent tutorial 
overview), as it was suggested as a feasible mechanism to ex¬ 
plain long-lived oscillations in the 2D spectra of several nat¬ 
ural light harvesting complexes and a photosynthetic reaction 
centeri^iiS. The requirement of exciton-vibrational resonance 
is readily satisfied in such systems, given their numerous exci- 
tonic bands and rich vibrational structures. Incoherent models 
based upon correlated fluctuations were not ruled out though. 
Our work provides a quantum mechanical foundation for en¬ 
hanced energy transfer based on vibronic coupling. As re¬ 
cently demonstrated, this mechanism is not limited to natural 
light harvesting, vibronic coupling is also of key importance 
in photovoltaic devicesi^. 

Methods 

Polarization controlled 2D electronic spectroscopy. In 

2D electronic spectroscopy, three ultrashort laser pulses gen¬ 
erate an optical response of a molecular ensemble, which 
is spectrally resolved along both absorption {a>i) and detec¬ 
tion (< 413 ) frequencies within the laser pulse spectrum. The 
absorption frequency a>i is obtained by precise scanning of 
the time delay between the first two pulses and subsequent 
Fourier transformation (fi —» a>i). In detection, the signal 
is spectrally dispersed, leading directly to the detection fre¬ 
quency W 3 . Varying time delay f 2 between pulses 2 and 3 
provides information about evolution of the system on a fem¬ 
tosecond timescale^—. In order to retrieve the purely ab¬ 
sorptive part, the signal induced by pulses 1-3 is detected in 
a heterodyned fashion by interfering it with a phase-stable 
local oscillator pulse (LO). Polarization control is achieved 
by the combination of A/4- wave plates and wire grid polar¬ 
izers for each of the laser beams to select the desired po¬ 
larization with high accuracy. Polarization-resolved 2D ex¬ 
periments change the relative contributions of distinct path¬ 
ways depending on the polarization of the laser pulses, orien¬ 
tation of the transition dipole moments and isotropy of the 
sample^. Rephasing spectra were acquired with a polar¬ 
ization sequence of (90°, 0°, 90°, 0°) for pulses (1,2,3, LO), 
in contrast to non-rephasing spectra, where the time order¬ 
ing of the first two pulses is reversed, leading to a polariza¬ 
tion sequence of (0°, 90°, 90°, 0°). The polarization scheme 
used for rephasing spectra (Fig. Id) shows 0° was defined 
to be parallel to the sample flow direction, depicted as a 
white arrow in Fig. la. For a macroscopically aligned sample, 
this particular polarization sequence selects pathways stem¬ 
ming from interband coherences and vibronic mixin g'^i*^ , dis- 
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cussed throughout the paper, while pathways with all-parallel 
transition dipole moments such as ground state bleach, stimu¬ 
lated emission, excited state absorption and also vibrational 
wave packet excitation are suppressed. For the details re¬ 
garding the experimental methods, see Supplementary Note 1. 
To subtract the non-oscillatory signals from 2D spectra, we 
employed a decay associated spectra analysis^, where the 
population decays were htted by a sum of three 2D-spectra 
with individual decay constants. The aJ 2 -maps in Fig. 2 were 
obtained using Fourier transformation (f 2 ^ 0 ) 2 ) with zero¬ 
padding up to 2^ data points. All measurements were carried 
out at room temperature. 

References 

*Van Amerongen, H., Valkunas, L. & Van Grondelle, R. Photosynthetic ex- 
citons (World Scientific, Singapore, 2000). 

^Blankenship, R. E. Molecular mechanisms of photosynthesis (Blackwell 
Science, Oxford, 2002). 

^Renger, T., May, V. & Kiihn, O. Ultrafast excitation energy transfer dy¬ 
namics in photosynthetic pigment-protein complexes. Physics Reports 343, 
137-254 (2001). 

"^Huelga, S. F. & Plenio, M. B. Vibrations, quanta and biology. Contemp. 
Phys. 54, 181-207 (2013). 

^Jonas, D. M. Two-dimensional femtosecond spectroscopy. Anna. Rev. Phys. 
Chem. 54 , 425^63 (2003). 

^Engel, G. S. et al. Evidence for wavelike energy transfer through quantum 
coherence in photosynthetic systems. Nature 446 , 782-786 (2007). 

^Dostal, J., Mancal, T., Vacha, F., Psencik, J. & Zigmantas, D. Unraveling 
the nature of coherent beatings in chlorosomes. J. Chem. Phys. 140 , 115103 
(2014). 

^Collini, E. etal. Coherently wired light-harvesting in photosynthetic marine 
algae at ambient temperature. Nature 463, 644-647 (2010). 

^Romero, E. et al. Quantum coherence in photosynthesis for efficient solar- 
energy conversion. Nature Phys. 10 , 676-682 (2014). 

^^Fuller, F. D. et al. Vibronic coherence in oxygenic photosynthesis. Nature 
Chem. 6, 706-711 (2014). 

^^Chin, A. W. et al. The role of non-equilibrium vibrational structures in 
electronic coherence and recoherence in pigment-protein complexes. Na¬ 
ture Phys. 9, 113-118 (2013). 

^^Plenio, M. B., Almeida, J. & Huelga, S. F. Origin of long-lived oscillations 
in 2D-spectra of a quantum vibronic model: electronic versus vibrational 
coherence. J. Chem. Phys. 139 , 235102 (2013). 

^^Chin, A. W., Huelga, S. F. & Plenio, M. B. Coherence and decoherence 
in biological system: principles of noise assisted transport and the origin of 
long-lived coherences. Phil. Trans. Act. Roy. Soc. A 370 , 3638-3657 (2012). 
^"^Kolli, A., O’Reilly, E. J., Scholes, G. D. & Olaya-Castro, A. The fundamen¬ 
tal role of quantized vibrations in coherent light harvesting by cryptophyte 
algae. J. Chem. Phys. 137 , 174109 (2012). 

^^Tiwari, V., Peters, W. K. & Jonas, D. M. Electronic resonance with anticor¬ 
related pigment vibrations drives photosynthetic energy transfer outside the 
adiabatic framework. PNAS 110 , 1203-1208 (2013). 

^^Lee, H., Cheng, Y.-C. & Fleming, G. R. Coherence dynamics in photosyn¬ 
thesis: protein protection of excitonic coherence. Science 316 , 1462-1465 
(2007). 

^^Ishizaki, A., Calhoun, T. R., Schlau-Cohen, G. S. & Fleming, G. R. Quan¬ 
tum coherence and its inteiplay with protein environments in photosyn¬ 
thetic electronic energy transfer. Phys. Chem. Chem. Phys. 12 , 7319-7337 
( 2010 ). 

^^Hayes, D., Griffin, G. B. & Engel, G. S. Engineering coherence among 
excited states in synthetic heterodimer systems. Science 340 , 1431-1434 
(2013). 

^^Christensson, N. et al. High frequency vibrational modulations in two- 
dimensional electronic spectra and their resemblance to electronic coher¬ 
ence signatures. J. Phys. Chem. B 115 , 5383-5391 (2011). 
^^Caycedo-Soler, F., Chin, A. W., Almeida, J., Huelga, S. F. & Plenio, M. B. 
The nature of the low energy band of the Fenna-Matthews-Olson complex: 
vibronic signatures. J. Chem. Phys. 136 , 155102 (2012). 


^^Christensson, N., Kauffmann, H. F, Pullerits, T. & Mancal, T. Origin of 
long-lived coherences in light-harvesting complexes. J. Phys. Chem. B 116 , 
7449-7454 (2012). 

^^Heijs, D.-J., Dijkstra, A. G. & Knoester, J. Ultrafast pump-probe spec¬ 
troscopy of linear molecular aggregates: effects of exciton coherence and 
thermal dephasing. Chem. Phys. 341, 230-239 (2007). 

^^Wiirthner, F, Kaiser, T. E. & Saha-Moller, C. R. J-aggregates: from 
serendipitous discovery to supramolecular engineering of functional dye 
materials. Chem. Int. Ed. 50, 3376-3410 (2011). 

^"^Eisele, D. M. et al. Robust excitons inhabit soft supramolecular nanotubes. 
PNAS 111 , E3367-E3375 (2014). 

^^Yuen-Zhou, J. et al. Coherent exciton dynamics in supramolecular light¬ 
harvesting nanotubes revealed by ultrafast quantum process tomography. 
ACS Nano 8, 5527-5534 (2014). 

^^Qiao, Y. et al. Nanotubular J-aggregates and quantum dots coupled for effi¬ 
cient resonance excitation energy transfer. ACS Nano 9 , 1552-1560 (2015). 

^^Von Berlepsch, H. & Bottcher, C. The morphologies of molecular cyanine 
dye aggregates as revealed by cryogenic transmission electron microscopy 
(World Scientific, Singapore, 2012). 

^^Von Berlepsch, H., Kirstein, S. & Bottcher, C. Effect of alcohols on J- 
aggregation of a carbocyanine dye. Langmuir 18 , 7699-7705 (2002). 

^^Von Berlepsch, H., Kirstein, S. & Bottcher, C. Controlling the helicity of 
tubular J-aggregates by chiral alcohols. J. Phys. Chem. B 107 , 9646-9654 
(2003). 

^^Von Berlepsch, H. et al. Supramolecular structures of J-aggregates of car¬ 
bocyanine dyes in solution. J. Phys. Chem. B 104 , 5255-5262 (2000). 

^^Von Berlepsch, H. et al. Stabilization of individual tubular J-aggregates by 
poly(vinyl alcohol). J. Phys. Chem. B 107 , 14176-14184 (2003). 

^^Eisele, D. M. et al. Utilizing redox-chemistry to elucidate the nature of 
exciton transitions in supramolecular dye nanotubes. Nature Chem. 4, 655- 
662 (2012). 

^^Milota, F. et al. Vibronic and vibrational coherences in two-dimensional 
electronic spectra of supramolecular J-aggregates. J. Phys. Chem. A 117 , 
6007-6014 (2013). 

^“^Aydin, M., Dede, O. & Akins, D. L. Density functional theory and 
Raman spectroscopy applied to structure and vibrational mode analysis 
of 1,1 ’,3,3’-tetraethyl-5,5’,6,6’-tetrachloro-benzimidazolocarbocyanine io¬ 
dide and its aggregate. J. Chem. Phys. 134 , 064325 (201 1 ). 

^^Rich, C. C. & McHale, J. L. Resonance Raman spectra of individual exci- 
tonically coupled chromophore aggregates. J. Phys. Chem. C 117 , 10856- 
10865 (2013). 

^^Butkus, V., Zigmantas, D., Abramavicius, D. & Valkunas, L. Distinctive 
chai'acter of electronic and vibrational coherences in disordered molecular 
aggregates. Chem. Phys. Lett. 587, 93-98 (2013). 

^^Womick, J. M. & Moran, A. M. Exciton coherence and energy transport 
in the light-harvesting dimers of allophycocyanin. J. Phys. Chem. B 113 , 
15747-15759 (2009). 

^^Womick, J. M. & Moran, A. M. Vibronic enhancement of exciton sizes 
and energy transport in photosynthetic complexes. J. Phys. Chem. B 115 , 
1347-1356 (2011). 

^^Del Rey, M., Chin, A. W., Huelga, S. F. & Plenio, M. B. Exploiting struc¬ 
tured environments for efficient energy transfer: the phonon antenna mech¬ 
anism. J. Phys. Chem. Lett. 4, 903-907 (2013). 

^^Hochstrasser, R. M. Two-dimensional IR-spectroscopy: polarization 
anisotropy effects. Chem. Phys. 266, 273-284 (2001). 

'^^Read, E. L. et al. Cross-peak-specific two-dimensional electronic spec¬ 
troscopy. PNAS 104 , 14203-14208 (2007). 

^^Mukamel, S. Principles of nonlinear optical spectroscopy (Oxford Univer¬ 
sity Press, Oxford, 1995). 

^^Perlfk, V. et al. Vibronic coupling explains the ultrafast carotenoid-to- 
bacteriochlorophyll energy transfer in natural and artificial light harvesters. 
J. Chem. Phys. 142, 212434 (2015). 

'^‘^Schrbter, M. et al. Exciton-vibrational coupling in the dynamics and spec¬ 
troscopy of Frenkel excitons in molecular aggregates. Phys. Rep. 567 , 1-78 
(2015). 

'^^Killoran, N., Huelga, S. F. & Plenio, M. B. Enhancing light-harvesting 
power with coherent vibrational interactions: a quantum heat engine pic¬ 
ture. Preprint at http://arxiv.org/abs/1412.4136 (2015). 

'^^Didraga, C. etal. Structure, spectroscopy, and microscopic model of tubular 
carbocyanine dye aggregates. J. Phys. Chem. B 108, 14976-14985 (2004). 


7 


^^Halpin, A. et al. Two-dimensional spectroscopy of a molecular dimer un¬ 
veils the effects of vibronic coupling on exciton coherences. Nature Chem. 
6, 196-201 (2014). 

"^^Chenu, A. & Scholes, G. D. Coherence in energy transfer and photosynthe¬ 
sis. Anna. Rev. Phys. Chem. 66, 69-96 (2015). 

^^Falke, S. M. et al. Coherent ultrafast charge transfer in an organic photo¬ 
voltaic blend. Science 344 , 1001—1005 (2014). 

^'^Brixner, T., Mancal, T., Stiopkin, I. V. & Fleming, G. R. Phase-stabilized 
two-dimensional electronic spectroscopy. J. Chem. Phys. 121 , 4221—4236 
(2004). 

^'Augulis, R. & Zigmantas, D. Two-dimensional electronic spectroscopy 
with double modulation lock-in detection: enhancement of sensitivity and 
noise resistance. Opt. Express 19 , 13126-13133 (2011). 

^^Augulis, R. & Zigmantas, D. Detector and dispersive delay calibration is¬ 
sues in broadband 2D electronic spech'oscopy. J. Opt. Soc. Am. B 30 , 1770- 
1774 (2013). 

Acknowledgements 

The authors would like to thank Valentyn I. Prokhorenko 
for help in 2D-DAS analysis. C.N.L. and J.H. acknowledge 
funding by the Austrian Science Fund (FWF): START project 
Y 631-N27 and by COST Action CM 1202 - PERSPECT- 
H20. J.L., E.C.-S., S.E.H. and M.B.P. acknowledge fund¬ 
ing by the EU STREP PAPETS and QUCHIP, the ERC 


Synergy Grant BioQ, the Deutsche Eorschungsgemeinschaft 
(DEG) within the SEB/TRR21 and an Alexander von Hum¬ 
boldt Professorship. J.P acknowledges funding by the Span¬ 
ish Ministerio de Economia y Competitividad under Project 
No. PIS2012-30625. D.P. and D.Z. acknowledge funding by 
the Swedish Research Council and Knut and Alice Wallenberg 
Eoundation. 

Author contributions 

D.P, D.Z. and J.H. designed and conducted experiments. 
H.v.B. was responsible for sample preparation, structural 
characterization and Raman measurements. J.L., E.C.-S., 
C.N.L., D.P, J.P, and J.H. analyzed the data. J.L., E.C.-S., 
S.E.H., J.H. and M.B.P. developed theory. All authors dis¬ 
cussed the results and wrote the manuscript. 

Competing financial interests 

The authors declare no competing financial interests. 

Correspondences and requests for materials should be ad¬ 
dressed to J.H. (juergen.hauer@tuwien.ac.at). 


Supplementary Information 

I. EXPERIMENT 
A. Sample preparation 

The monomer, tetrachlorobenzimidacarbocyanine chromophore with two attached hydrophobic octyl groups (EEW- 
Chemicals, Wolfen, Germany) was dissolved in 10“^ M NaOH solution to achieve a concentration of 10“^ M. The solution 
was then stirred in the dark for several hours. Subsequently, polyvinyl alcohol (PVA) of molecular weight ~ 130000 was added 
in 1;10 w/w ratio (monomer;PVA) to slow down the formation of aggregate bundles during the storage of dye solutions. More¬ 
over, the adsorbed PVA chains^ obviously prevent the reassembly of double-layered into single-layered tubes upon bundling. 
This effect was observed recently for another derivative^ (C8S3) of the present dye. The individual tubular aggregates degrade 
in that case into single-layered tubes, which is accompanied by a dramatic change of absorption spectra. Similar effects were 
not observed for C803 when PVA is present. In particular, the aggregate solutions prepared in the described way were stable for 
approximately ten days when stirred continuously. Without stirring the spectral signature of the double-layered tubes retained 
even after 12 weeks of storage^. Eor 2D experiments, we additionally diluted the sample with 10“^ M NaOH to obtain optical 
density below 0.3 at 598 nm at a path length of 200 /im. 

A total sample volume of approximately 10 ml was circulated through the U-shaped wire-guided jeti by a peristaltic pump 
(Masterflex C/L) with a flow speed optimized for film stability. Solvent evaporated from the recollecting container was refilled 
every 4 hours during the course of 13 hour measurement. 


B. Data acquisition 

Passively stabilized 2D spectroscopy was described in detail elsewhere^. Briefly, a home-built non-collinear optical parametric 
amplifier (NOPA) seeded by 180 fs pulses at 1030 nm from PHAROS (Light Conversion Ltd) was tuned to generate ~ 16 fs pulses 
(80 nm full width at half maximum) centered at 580 nm. The NOPA output was split into four pulses and arranged in the so-called 
boxcar geometry. Waiting time t 2 was controlled by a mechanical translation stage (PI), whereas coherence time t\ was scanned 
by inserting a pair of fused silica wedges into the first two pulses. All four pulses were focused and overlapped in the sample. 
The first three generated a third order nonlinear optical response which is emitted in the photon echo phase-matched direction. 
This signal was heterodyned with an attenuated fourth pulse, called local oscillator (LO). The resulting interference pattern was 
spectrally resolved and detected by a CCD camera (PIXIS, Princeton Instruments). Most of the scatter was eliminated by the 
double-frequency lock-in modulation of the first two pulses^. The polarization of each pulse was controlled by the combination 
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-o-PP 
— proj. 2D 


16.0 16.5 17.0 17.5 18.0 18.5 

Detection frequency cOg /10^ cm'^ 


FIG. 4. Pump-probe and projected 2D signal. Projection of phased polarization-controlled 2D spectra (blue) to all-parallel pump-probe 
(light blue) at t 2 = 180 fs. 


of AjA wave plates and wire grid polarizers (contrast ratio > 800). The accuracy of the polarization angle was estimated to be 
+ 1 °, where the unwanted signals were typically suppressed by a factor of ~80 for the selected polarization sequence. 

To prevent degradation of the sample, the power and repetition rate of the laser were set to 200pJ/pulse and 40 kHz, respec¬ 
tively. Spectral resolution of ~35 cm ' for the detection frequency W 3 was determined by the grating, the number of CCD pixels 
and Fourier filtering of the signal during the standard analysis procedure. Coherence time was scanned from -300 fs to 384 fs 
in 1.5 fs steps, providing ~43 cm“' spectral resolution of absorption frequency a>i. Waiting time steps of 12 fs were sufficient to 
resolve oscillatory features up to 1350 cm“* with ~35 cm“' resolution along a> 2 - 


C. Polarization-controlled 2D-ES 

The strength of 2D signals is determined by the scalar products of molecular transition dipole moments and pulse polarizations. 
To take advantage of i) the preferential orientation of the J-aggregate (from here on referred to as C803) along the flow direction 
of the jet and ii) mutually perpendicular transition dipole moments of bands 1(2) and 3 of C803, we designed a polarization 
scheme selective for interband coherences. This is similar to the case of an isotropic sample discussed both theoretically^ and 
experimentally2i2. In the presented experiments, the polarization scheme for rephasing signals reads (90,0,90,0) for beams 
1-4, respectively. The first and third pulses, polarized orthogonal (90) to the jet’s flow direction, interact with bands 1-3. The 
second pulse, polarized parallel ( 0 ) to the jet’s flow direction, interacts preferentially with bands 1 and 2 , due to the negligible 
transition dipole moment of band 3 along this direction. The polarization scheme for non-rephasing spectra reads (0,90,90,0), 
as the ordering of the first two pulses is reversed. These polarization schemes restrict oscillatory signals induced by interband 
coherences to the lower cross peak in rephasing spectra (R31) and the lower diagonal peak in non-rephasing spectra (Nil), as 
shown in Figures 2 and 3 of the main text. Non-oscillatory 2D signals were subtracted prior to Fourier transformation t 2 coi- 

The polarization-controlled 2D spectra were phased to pump-probe data where pump and probe pulses were polarized in 
parallel. This procedure is not rigorously correct because in polarization-controlled 2D-ES the first two pulses have different 
polarization directions while in pump-probe the first two interactions derive from the same pump pulse which naturally means 
parallel interactions. In other words, the projection slice theorem is strictly speaking not valid for the experiments presented 
her&iS. Despite this discrepancy, one can still satisfactorily phase polarization-controlled 2D spectra to all-parallel pump-probe 
as shown in Supplementary Figure |4] One explanation of this is leakage of the much stronger all-parallel signals through the 
crossed polarizers, meaning that the all-parallel signal still dominates the non-oscillatory part of the (90,0,90,0) 2D-signal. In 
this work, we decided to phase polarization-controlled 2D data to parallel pump-probe data. We note that the imperfection in 
phasing parameters only affects the lineshapes of the real and imaginary part of 0)2 maps, but preserves their amplitude-maps in 
both lineshape and magnitude. Hence, the difficulties in phasing polarization-controlled 2D spectra discussed above do not affect 
the conclusions drawn in the main part of the main text, which were based on 0)2 amplitude-maps. To this end, we found that 
arbitrary and large changes of the phasing parameters do not alter 012 amplitude-maps shown in Figure 2 of the main text (results 
not presented). It is noted that sophisticated phasing techniques based on heterodyned transient grating instead of pump-probe 
offer a correct method to phase crossed-polarization 2D signalsil. 
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II. THEORY 

A. A vibronic model for bands 1 and 3 of C803 

In the following the vibronic model used to describe bands 1 and 3 of C803 and simulate 2D spectra is described. We consider 
coherent interaction of bands 1 and 3 with the intramolecular vibrational modes of frequency tLV\ 668 cm ', which is quasi¬ 
resonant with the exciton energy splitting between bands 1 and 3. The environmental noise induced by background phonons (a 
phonon bath) is modeled by a Markovian quantum master equation. 


f. Hamiltonian 

The electronic Hamiltonian of C803 that consists of a network of cyanine dye molecules is described by 

h, = 2 \^a} (^orl "I" ^ ^ l^or) (1) 

a 

( 2 ) 

k 

where |e„) represents the excited state of site a (or molecule a), Ea denotes the site energy including electronic and reorganization 
energies, and Jap the electronic coupling between sites a and p. The diagonalization of the electronic Hamiltonian Hg gives rise 
to the exciton states \k) - \ea) {ea\k) associated with the exciton energies Qj., where bands 1 and 3 are denoted by |1) and |3), 
respectively: {ea\k) 6 R for Ea, Jap 6 R- 

The vibrational modes with frequency hvi » 668 cm ' are described by a set of harmonic oscillators 

Hy ^ 'ynvialag, (3) 

a 

where and Ug represent the creation and annihilation operators, respectively, of the intramolecular vibrational mode of site a. 
The interaction between vibrations and the electronic excitation of molecules is modeled by 

He-v = ^ f^vi V^T k«) <eQ.| (al + Oa), (4) 

a 

where si denotes the Huang-Rhys factor of the vibrational modes. In the exciton basis {|k)), the interaction Hamiltonian is 
represented by 

Hg-v — hv\ \k) {l\'^{k\ea){ea\l) {al + Qa), (5) 

k,l a 

where the diagonal terms {k - 1) lead to adiabatic surfaces in the electronic excited states, called vibrons, while the non-diagonal 
terms {k + 1) induce coherent transition between different excitons mediated by exciton-vibrational couplings. 

In this work, we are interested in the coherent interaction of bands 1 and 3 with the quasi-resonant vibrational modes of 
frequency vi, which is described by the following cross term Hg-y in Eq. (|5]l 

Hg.y = hvx V^(|1><3| + |3)<i|)(a[ + ai), (6) 

where ai - N 2a (lka)(ea|3)aa describes an effective vibrational mode with frequency vi. Here is introduced to normalize the 
effective vibrational mode, such that [5i, aj] = 1, leading to an effective Huang-Rhys factor S \ - sijN'^. This implies that for a 
given Huang-Rhys factor ii, the effective Huang-Rhys factor 5 1 is increased as the spatial overlap {l|ea)(ea|3) between excitonic 
wavefunctions of bands 1 and 3 increases, leading to smaller N and larger S i = sxjN'^. The effective Hamiltonian of bands 1 
and 3 coupled to the effective vibrational mode is then described hy H - Hg + Hy + Hg-y, where Hg - HQ.i |1) {1| -i- |3) {3| 

and Hy - Jivia\ai. We note that the vibrational energy hvi » 668 cm“* is higher than the thermal energy kgT » 208 cm ' at 
room temperature T - 300 K, implying that the thermal state of the vibrational mode is well approximated by its ground state. 
In addition, when exciton-vibrational couplings are sufficiently small, the light-induced vibrational excitation of overtones is 
negligible due to the small Franck-Condon factors. This is the case for C803, where Nil and R31 in 2D spectra can be well 
described within a subspace spanned by {|go), ki), |lo), |li), |3o))- Here, |ko) and |ki) denote the vibrational ground and first 
excited states of an electronic state \k), respectively, i.e. (Hg+Hy) |k;) = fi{Q.k + lvi) |k;), where |g) represents the electronic ground 
state with = 0. In this scenario, {|lo), |3o)) can be directly excited by light from the ground state |go), while {|gi), |li)) has 
an extremely low transition probability due to small Franck-Condon factors. Nonetheless, {|gi), |li)) can be populated through 
exciton-vibrational coupling vi leading to transition from |3o) to 11 1 ), and subsequently to |gi) via emission. The coherent 
transition between 1 3o) and 111 ) requires resonance between vibrational frequency vi and exciton energy splitting AQ 31 = D 3 - D 1 
between bands 1 and 3, i.e. AD 31 a; vi. 
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2. Decoherence 


In addition to the coherent interaction of bands 1 and 3 with the effective vibrational mode 5i, we consider electronic decoher¬ 
ence induced by background phonons. We characterize the decoherence by two dynamical processes, i) the incoherent population 
transfer between excitons, called exciton relaxation, and ii) the pure dephasing noise that destroys electronic coherence without 
exciton population transfer. In addition, we consider iii) relaxation of the effective vibrational mode. 

We assume that each cyanine dye molecule is coupled to an independent phonon bath. The Hamiltonian of the background 
phonons is given by //ph = hv^b''^b^ with the interaction Hamiltonian \ca}{ea\ + b^) between molecules 

and phonons, where and b^ denote the creation and annihilation operators, respectively, of a background phonon mode f. 
Here ga^ represents the exciton-phonon coupling between site a and phonon mode which satishes ga^gp^ - 0 for all P + a, 
implying that when site a is coupled to the phonon mode ^ with gaf + 0, all the other sites p are decoupled from the mode 
with gp^ - 0. For the sake of simplicity, we assume that there is no degeneracy in the exciton energies Q.k, which leads to 
a relatively simple form of a Markovian quantum master equation. This condition is satished even if the exciton energies are 
close to degeneracy unless they are strictly degenerate, which is satished for bands 1 and 3 of our interest. The inhuence of the 
background phonons on the vibronic system consisting of bands 1 and 3 with the effective vibrational mode is then described by 
a Markovian quantum master equatioi>i^ 

^P(0 = + DMt)] + r)a[p{t)] + DMt)l (7) 

at n 

where p(t) denotes the reduced vibronic state, while Dr[p{t)], Dd[p{t)] and Dv[p(t)] describe exciton relaxation, pure dephasing 
noise and relaxation of the effective vibrational mode, respectively. 


i. Exciton relaxation 

Here 2)r[p(f)] describes exciton relaxation 

n[p(f)] = 2 Z lAa(aj)p(t)Alioj) - i{A^(cu)A„(cu),p(f))), (8) 

W#0 Q- ' ' 

with Aflj.; = Q.k - Q.I denoting the exciton energy splitting between \k} and \l), Aa{co) - Ailki) {l\ea) {ea\k)\l) {k\ for 

0 ) + 0, leading to incoherent transition from \k) to |f), where 6 {i, j) denotes the Kronecker delta dehned by 6 (i, j) = 1 if f = j and 
d{i, j) - 0 otherwise. In Eq. ®, jaaitA) is dehned by 


yaaioj) = 27Tjaia>)in(L0) + 1), (9) 

with n{a)) - {exp(fia)IksT) - 1)“* representing the Bose-Einstein distribution function at temperature T, J'aioS) is the spectral 
density of site a dehned by J'aitA) = - u^) if w > 0 and J'aiop - -,!7a(-w) otherwise. Here (5(v) represents the Dirac 

X OO 

^ dx6(x) — 1. 

ii. Pure dephasing noise 

^d\p(ty\ in Eq. (|7]i describes the pure dephasing noise 

2)^[p(0] = Z (A.(0)p(f)Ai(0) - ^{Ai(0)A„(0),p(f))j, (10) 


where Aq.(0) = Yjk K^ko-)!^ k) (^1 destroys electronic coherence without changing exciton populations dehned by {Tr[(k| p(f) |k)]). 

By substituting electronic coherences |g)(l|, |g){3| and |1)(3| to the dissipators 2)r[p(f)] and 2)rf[p(f)] in Eqs. (|8j and ( fTOl i. 
one can obtain the following electronic decoherence rates E^i, rg 3 and ri 3 of the coherences |g) (1|, |g) {3| and |1) {3| 


r«i 

r«3 


ri3 



+ r«i> 

m 



+ r«3, 

m 



-3^ 


/#i m 


( 11 ) 

( 12 ) 


(13) 
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where denotes the incoherent population transfer rate from band k to I 

jk^i = jaaiASlki) K/|eQ-)<eQ.|^)|^ > 0, (14) 

a 

while -ygi and 7^3 represent the pure dephasing rates of the coherences |g) (1| and |g) (3|, respectively, 

a 

a 

and 713 represents the pure dephasing rate of the inter-exciton coherence |1) {3| between bands 1 and 3 

a 

These results imply that the inter-exciton dephasing rate 713 should be lower than the sum of the other dephasing rates 7^1 and 
7 g 3 when there is a spatial overlap between excitonic wavefunctions of bands 1 and 3 

ri3 =r«i +7^3 - ^Klkff)|^7aQ.(0)K3|ea)|^ <7i;i H-7g3, (18) 


(15) 

(16) 


with jaaiO) > 0 for all a, the equality 713 = 7^1 -t- 7^3 holds if and only if there is no spatial overlap between excitonic 
wavefunctions, i.e. Kl|ea)P K 3 |ea)|^ = 0 for all a, or the spectral densities Sfaio^) of the molecules shared by bands 1 and 3 do 
not induce pure dephasing noise by yaa(O) - 0 for all sites a satisfying |(l|eQ.)|^ K 3 |ea)|^ 4^ 0 . This implies that even if each 
molecule is coupled to an independent phonon bath, the spatial overlap between excitonic wavefunctions can reduce the inter- 
exciton dephasing rate 713. Here the independent phonon baths of the molecules shared by excitons effectively form a common 
phonon bath coupled to both excitons, leading to a partial dephasing-free subspace. For instance, if bands 1 and 3 have perfect 
spatial overlap, i.e. Kl|ea)|^ = K 3 |ea)P for all a, while the orthogonality between them is satisfied by the phases of (Ika) {ea\'i), 
i.e. ( 1 | 3 ) = 2 q. {lka)(eQ.| 3 ) = 0 , the inter-exciton dephasing rate 713 will become zero, as each Aq.( 0 ) = Yjk K^ka)l^ k) (^1 = 
Klko-)!^ (| 1 ) ( 1 | + | 3 ) ( 31 ) + K^ka)P k) (^1 forms a dephasing-free subspace | 1 ) ( 1 | + | 3 ) ( 3 | of bands 1 and 3 . Since band 1 
is localized on the inner layer of C 803 , while band 3 is delocalized on both the inner and outer layersi^, there is a partial spatial 
overlap between excitonic wavefunctions, leading to 713 < 7^1 + 7^3. The spatial overlap is also required for a non-zero value of 
the effective Huang-Rhys factor S 1, which is responsible for the long-lived beating signals observed in the experiment, as will 
be discussed later. 

In addition, the inter-exciton dephasing rate 713 has a non-zero lower bound when the dephasing rates jgi and 7^3 are different 
in magnitude. The dephasing rates in Eqs. (fTSt-tflTTi can be expressed as 7^1 = |vip, 7^3 = IV3P and 713 = |vi - V3P with the 
real vectors i4 defined by i4 = where Wk is a real vector with elements \{k\ea)?' > 0 representing the delocalization 

of an exciton state k) in the site basis (ka)), while 7 is a diagonalized matrix with elements 7aa(0) > 0, leading to a positive 
matrix 7'^^ defined by 7 = pj-om jpe triangle inequality, |vi - V3I + IV3I > |vi| and |vi - V3I -H |vi| > |v3|, the inter-exciton 

dephasing rate 713 is bounded from below by ( ^Jjgi - ^Jjgi )^, leading to ( yjjgx - yjjgi < 713 < Jgi + Jg3- Therefore, the 
electronic decoherence rate Tn of the inter-exciton coherence | 1 ) { 3 | is constrained by 


Ml M3 




y/2 

/ 




y/2 

•; 


2 


< ri3 < Fgi H- rg3. 


(19) 


with jgk — Fg^. - i Yii+k 7k^i from Eqs. (fTTT l and (fT^ . Here the population transfer rates ( 71 ^; and 73 ^/) and electronic 
decoherence rates (Fgi and Fg 3 ) can be estimated using experimentally measured 2D spectra, which will be discussed later. 


iii. Relaxation of quasi-resonant vibrations 

Einally, fDi,[p(f)] in Eq. (|7|i describes the relaxation of the effective vibrational mode 

®v[p(0] = 7 v{n(v\) + 1) (2aip(f)a| - {a|ai,p(f))) + 7vn{v\} (2a\p{f}ai - {aia\,p{f}}) . (20) 

Since n(vi) =» 0.04 at room temperature T - 300 K due to the high vibrational energy hvi » kgT, Eq. (l20l) can be reduced to 

D,.ip(t)] » 7„(2aip(f)d| - {a|ai,p(f))), (2i) 

which describes the dissipation of the vibrational mode with the rate of 7 ,,. 
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FIG. 5. Feynman diagrams contribnting to the beating signals in Nil and R31 represented in uncoupled state basis, a. The stimulated 
emission diagram contributing to the beating signals in N1 1 . Here time runs upwards and the electronic transitions induced by light are denoted 
by arrows: 0 and 90 denote the polarization of light, parallel and normal to the longitudinal axis of C803, respectively (cf. Figure 1 in the main 
text). The time interval between pulses is called coherence time f,, waiting time t 2 , and rephasing time for the first and second, the second 
and third, and the third excitation pulse and the emerging signal, respectively. The Fourier transform along fj and leads to the absorption 
and detection frequencies denoted by Wi and respectively, b-e. The stimulated emission and ground state bleaching diagrams contributing 
to the beating signals in R31. In a-d, grey shaded waiting time periods during t 2 highlight vibronic coherences in the electronic excited states. 
In e, on the other hand, the vibronic system is in the electronic ground state during t 2 . 


3. The response function for N11 

Here we derive the response function for the beating signals in Nil, which is a diagonal peak in non-rephasing spectra 
centered at (wi, W 3 ) « 

In Fig. |5^, the Feynman diagram contributing to the beating signals in N11 after the employed (0,90,90,0) excitation is 
displayed. As the thermal state of the effective vibrational mode at room temperature is well approximated by its ground state 
{hvi » kgT), the initial state of the vibronic system is given by |go)(^ol- After excitation to |lo){gol by the first pulse, the 
dynamics of |lo) {gol during coherence time fi is governed by a time evolution super-operatordetermined by the quantum 
master equation in Eq. © 




( 22 ) 


for which the Fourier transform is given by 





1 

i{coi - Qi) - Fgi 


HoX^ol, 


(23) 


where the prefactor -(iicoi - Qi) - F^i) ' determines the lineshape of Nil along the wi-axis, which is centered at u)\ - Qi 
with a linewidth of 2rgi. By the second pulse, | Iq) (gol becomes | Iq) (3o|, which evolves during waiting time f 2 into a mixture of 
|lo) (3o| and |lo) (lil, mediated by exciton-vibrational coupling, scaling with vi The time evolution of |lo) (3o| is formally 
expressed as 


miMW <3ol = Ho) <3ol = /(f 2 ) Ho) <3ol + g(f 2 ) Ho) dil, 
where is a non-Hermitian operator describing both the Hamiltonian dynamics and decoherence 


(24) 


K = (/An3i - Tb) |3o> <3o| + (ivi - Tv) Hi) (lil + ;Vi V^(|3o) <li| + H 1 ) <3ol). 


(25) 
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Here we evaluate /(f 2 ) in Eq. (l24l) . which describes the case that |lo) (3o| becomes |lo) {^ol by the third pulse, as shown in Fig.|5t. 
By diagonalizing the non-Hermitian operator K, one can show that /(f 2 ) is given by 


fiti) 

k=l ^ 


i + (-ir 


X — y 


- y)^ + 


exp 


^{x + y + (- 1 )* yjix-y)^ + 4z^^ t2 


where x - /AQ 31 - ri 3 , y = ivi - jv and z = ivi V^- Finally, |lo) {gol evolves during rephasing time t^, 


(26) 


(27) 


for which the Fourier transform leads to the lineshape-(!( w 3 - Qi) - F^i) ' ofNll along the a) 3 -axis. Therefore, the response 
function for N 11 is given by 


f^l^(Wi,f2,W3) - yipl4n~ 


1 


1 


Hcoi - ill) - r^i i(co 2 - Qi) - r,i 


m), 


(28) 


where ^ip denotes the transition dipole moment of band 1 for light polarized parallel to the longitudinal axis of C803, while 
fi 3 „ represents the transition dipole moment of band 3 for light polarized normal to the axis. This is due to the (0,90,90,0) 
polarization scheme employed for measuring non-rephasing spectra in the experiment, as schematically shown in Fig. |5^. It 
is notable that all the Feynman diagrams in Figs. |5^-e can be induced by (90,90,90,90) excitation where all the pulses are 
polarized normal to the longitudinal axis: band 1 can be excited or de-excited by both 0 and 90 polarizations, although with 
higher efficiency for light polarized at 0. For (90,90,90,90) excitation, the overall dipole strength in Eq. (|28]) is decreased 
to with as band 1 is mainly polarized along the longitudinal axis of C803, as shown in the linear dichroism 

spectrum in Figure 1 of the main text. This implies that the (0,90,90,0) polarization scheme for non-rephasing spectra enhances 
the signal-to-noise ratio when compared to the (90,90,90,90) excitation. Similarly, the signal-to-noise ratio of rephasing spectra 
is enhanced by (90,0,90,0) excitation. 

The lineshape function {i{a>i - Qi) - rgi)“^(/(w 3 - Qi) - F^i)"^ in Eq. (l28l l shows that N11 is centered at (wi, W 3 ) = (fli, fli) 
with a symmetric linewidth 2Tgi along wi- and W 3 -axes. When (a;i,cu 3 ) = (Qi,Qi), the lineshape function is reduced to F^^ 
implying that the amplitude of the Nil peak is proportional to F^^, which is decreased as the linewidth 2rgi increases. The time- 
dependent term /(f 2 ) in Eq. (l28l) describes the evolution of N11 during waiting time f 2 . In the absence of the exciton-vibrational 
coupling {S 1 = 0 ), /(f 2 ) is reduced to 

/(f2)l5.=o = (29) 


implying that the coherence |lo)(3ol oscillates with the frequency of the exciton energy splitting Afl 3 i and decays with the 
electronic decoherence rate ri 3 . Conversely, in the presence of the exciton-vibrational coupling {S 1 > 0), /(f 2 ) is expressed as 


1 

/ A 

1 , 

[/(An3i+5a>)-ri3+^y]f2 , ^ 

/ \ 

2 

1 ^J(x-y)^+Az^) 

2 

1 y/(x-y)^+4z^j 


li(vi- 6 oi)-y^,-Sy]t 2 


(30) 


where iSo) + 6 y - 2 “' [ yj(x - y)^ + Az^ - (x - y)], which satisfies 6 cj > 0 and dy > 0 for AQ 31 > vi and F^ > jv, which is the 
case for C803. There are several notable features that result from the vibronic coupling evident in Eq. (l30l l. i) The first term, 
proportional to e['(Af! 3 i+i 5 £<j)-ri 3 + 5 y]» 2 ^ oscillates with a frequency of Afl^j = Af 23 i - 1 - 6 a>, which is higher than the exciton energy 
splitting AQ 31 , and decays with the rate of ri 3 - 6 y, which is lower than the electronic decoherence rate F^ shown in Eq. ( |29] ). 
These are the characteristics of the vibronic coherence |lo)(3o|, where (3o| is one of the left eigenstates of K in the form of 
(3o| oc {3o| -I- ^(lil with 1^1 < 1. The vibronic eigenstate (3o| has a higher energy-level than {3o| due to the exciton-vibrational 
coupling, leading to AQ^j > AQ 31 (see Figure 3a in the main text). Additionally, the amplitude of |lo)(3o| in |lo){lil denoted 
by ^ leads to a longer lifetime than the coherence |lo){3ol that has no vibrational character, or in other words, the lifetime 
borrowing effect, ii) Conversely, the second term in Eq. ( l30t . proportional to e^A''i-S‘^'>-yv-Sy]t 2 ^ exhibits characteristics of the 
other vibronic coherence |lo) (Iil, where (ii| cx (li| - ^(3o| is the other left eigenstate of K. The second term oscillates with 
frequency Vj - vi - 6 a>, which is lower than the vibrational frequency vi due to the exciton-vibrational coupling (see Figure 3a 
in the main text). It also decays with the rate of y,, -H 6 y, which is higher than the vibrational decoherence rate yi. of |lo) (lil due 
to the amplitude of |lo) (lil in |lo) {3o| denoted by iii) We add that the vibronic states (3o| (3o| (lil and (Ii| oc {li| -^(3o| 

are the eigenstates of the non-Hermitian operator K in Eq. ( |25]) describing both Hamiltonian dynamics and decoherence, where 
^ depends on the parameters of the Hamiltonian as well as decoherence rates. These states are different from the eigenstates of 
the Hamiltonian H, which do not depend on decoherence rates, and their difference becomes non-negligible when the electronic 
decoherence rate ri 3 is comparable to or larger than the exciton-vibrational coupling vi as is the case for C803. 

By fitting experimental 2D spectra to the theoretical prediction of N11 and R31, which will be discussed later, we found that 
MQ 31 720cm-', hvi ^ 668 cm-', HTgi ^ 65 cm-', HTg^ ^ 150cm-', /jF^ ^ 80cm-', 5 1 = 0.0006 (cf. hvi VST ~ 16 cm-') 
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- Experiment 

- - Theory 
--Band 1 

- - Band 2 

- - Band 3 

- - Band 4 

- - Band 5 


FIG. 6. Absorption spectrum of C803. a, Absorption spectrum with light polarized parallel to the longitudinal axis of C803. Experimental 
and theoretical results are shown as a black solid line and a black dashed line, respectively. Theoretical results were modeled by a sum of 
Lorentzian functions, which describe bands 1-5 of C803. Each Lorentzian function is shown as a colored dashed line, b. Absorption spectrum 
with light polarized normal to the longitudinal axis of C803. Note that the vertical scales in a and b are different. 


and yv < (1 ps) '. The estimated electronic decoherence rates r,i and r,3 reproduce well the absorption spectrum of C803, as 
shown in Fig. | 6 ] where experimental and theoretical results are shown as a black solid line and a black dashed line, respectively. 
The theoretical results were modeled by a sum of the Lorentzian functions with linewidths 2Y'gk for k e {1,2,3,4,5), each of 
which describes the absorption of band k: each Lorentzian function is shown as a colored dashed line. The estimated values of the 
parameters lead to fidco ^ 1.6 cm“^ and My ^ 2A cm \ which are smaller than the experimental resolution of ~40cm '. This 
implies that for the case of C803, we can approximate and Vj by Af 23 i and vi, respectively, with 6a> ^ 0 and 6y ^ 0. More 
specifically, when the exciton-vibrational coupling is sufficiently small, such that vi V^T < l*Avi - Tnl with Avi = AQ 31 - vi, 
and the dissipation rate of the vibrational mode is negligible within the timescale of the total measurement time, i.e. y,, a; 0 , the 
response function determining N11 in Eq. (|28] ) is reduced to 




e 


(iAn3i-ri3)/2 






(31) 


with 62 representing the degree of vibronic mixing during waiting time t 2 


e 2 = iviV^(iAvi-ri3)-', (32) 

where the vibronic eigenstates {3ol and (Ii| are approximated by {3o| oc {3ol + 63 {li| and {li| oc (Ijl - 63 (3o|, respectively, with 
kaP 1 (in the main text, 63 was denoted by e for the sake of simplicity). It can be seen in Eq. (l32l) that |e 2 | increases as 
the exciton-vibrational coupling vi ■\/S~i increases or the detuning |Avi| = IAQ 31 - vi \ between exciton splitting and vibrational 
frequency decreases. This implies that the vibronic mixing of the coherences |lo) (3o| and |lo) (1 1 | requires resonance between 
excitons and vibrations and induces the observed long-lived beating signal in N11. In this respect, when |e 2 | decreases as a result 
of a high electronic decoherence rate Tn, the coherence |lo) {3o| generated by the second pulse (see Eig.|^) will decohere too 
quickly and thereby suppressing the vibronic mixing of |lo) (3o| and |lo) (1 1 | during waiting time 13 , which in turn will suppress 
the long-lived beating signal in N11. This is related to the fact that 63 is proportional to the exciton-vibrational coupling vi ViST 
and the amplitude of the long-lived component in Eq. (l3Tl l is proportional to e|. As such, when vi V^T < liAvi - EbI, 
the response function for N11 can be effectively described by two transitions between |lo) {3o| and |lo) (lil during waiting time 
f 2 , mediated by exciton-vibrational coupling vi i.e. |lo) {3o| ^ |lo) (111 —> ko) {3o|, within the timescale of the electronic 
decoherence rate E 13 , as shown in Eig.|5t. When the condition of vj Vs7 < |iAvi - r^l is not satisfied, the response function for 
N11 is represented by f 2 , W 3 ) = W 3 )(ivi V^)^" with the higher order terms proportional to (ivi 

which describe multiple transitions between |lo) (3ol and |lo) (1 1 | during 13 . 

In summary, when vi < I^Avi - T^l and y^ 0, the response function for N11 at (wi, 013 ) = (Qi, Qi) is given by 

Rig{t2) - [g(.Aa,.-r.3k3 + ^ (33) 

with 62 defined in Eq. (l32l l. The lineshape of Nil is symmetric along a>i- and W 3 -axes with a linewidth of 2rgi. These results 
are in line with the experimental observations shown in Eigures 2 and 3 of the main text. 


4. The response function for R31 

Here we provide the response function for the beating signals in R31, which is the cross peak in the rephasing spectra centered 
at (wi, W 3 ) =» (Q 3 , Qi). The response function for R31 can be derived using the same approach described above for N11. Here 
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we provide the results without derivation. Figs.|5})-e show the Feynman diagrams contributing to the beating signals in R31. In 
Figs-l^t-d, the vibronic system is in the electronic excited states during f 2 , while in Fig.|5^, the system is in the electronic ground 
state, each of which is called the stimulated emission (SE) and ground state bleaching (GSB) diagram, respectively. 

When vi ^/S~i < |/Avi - F 13 I, vi < |/Avi + Fgi - F^a] and jy « 0, which are satisfied for the case of C803, the contribution 
of the SE diagrams to R31 is approximated by 


R2g{0J\,t2, W 3 ) ~ t4pR3n 


1 


1 


-i(a>i - 0 . 3 ) - Fg 3 i{a )3 - Qi) - Fgi 
1 1 


e 


,(/An3i-ri3)f2 




(34) 


— ^^3) “ rg3 — Q3 + Avi) — Fgi I /(<x >3 — Oi) “ Fgi 


1 


. ^_g(<An3i-ri3)f2 


+ - 


1 


1 


1 


—i{coi — ^^ 3 ) “ r ^3 —i{coi — n 3 + Avi) — r^i j /(ct >3 “ “ Tgi 

with 61 representing the degree of vibronic mixing during coherence time ti 

ei = /vi V^(;Avi + Fgi - Fg3)^', 


(/An3i-ri3)f2^2 

^1 


(35) 


where AQ^ j and v' are approximated by AQ 31 and vi, respectively. More specifically, ei is associated with the transition between 
|go){3ol and |go)(lil during fi, while 62 is associated with the transition between |lo)(3o| and |lo){lil during f 2 . In Eq. (l34l l. 
the first term proportional to describes the transition |lo) (3ol —> |lo) (111 |lo) (3ol during t 2 (see Eig.|2)), the second term 
proportional to €162 describes the transition |go) (3ol ^ l^o)(lil during fi and the subsequent transition |lo)(lil ^ |lo)(3ol 
during t 2 (see Eig.|5};), and the last term proportional to describes the transition |go) (3ol —> Igo) (111 —> l.go) (3ol during fi (see 
Eig. |5}l). In the second and last terms, the lineshape function along the wi-axis contains (-i(wi - Q 3 + Avi) - F^i) ', which 
describes the presence of a sub-peak centered at co\ - - Avi = fli - 1 - vi < Q 3 with a linewidth of 2 Fgi, which is induced by 

exciton-vibrational coupling. However, due to the condition of |ei p « 1 and |e 2 p « 1, the first term in Eq. (l34l i determines the 
overall lineshape of R31, which is given by (-/(wi - ^ 3 ) -Fg 3 )“'(/(w 3 - fli) -Fgi)^' that is centered at (wi, CU 3 ) = (fl 3 , Qi) with 
the asymmetric linewidths of 2 Fg 3 and 2 F^i along wi- and W 3 -axes, respectively. 

The contribution of the GSB diagram to R31, with ground state coherence^ during f 2 , is given by 


R3g(Wi,f2,W3) : 


1 


1 


—i(toi — Q3) — Fg3 —i{toi — Q3 -I- Avi) — Tgi I \i{t 03 — Qi — Avi) — Fg3 i(w3 — fli) ~ Fgi 


1 


1 




(36) 


with 63 representing the vibronic mixing during 


€3 - -ivi -v^(-*Avi -i-Fji -Fg3) ’, ( 37 ) 

which is associated with the transition | 3 o) (gil ^ |li) (gil during shown in Eig.j^. Here the vibrational frequency vi in 
stems from the vibrational coherence |go)(^il in the electronic ground-state manifold and not the result of the approximation 
6aj^0. 

In summary, when vi -yST < |iAvi - F13I, vi V^T < |iAvi -1- Fgi - Fj3| and y,, « 0 , the response function forR 31 at (cui,a»3) = 
(^3, Gi) is given by 


R 2 g{t 2 ) + R3g(t2) ^ + e'''>'^e^(/ 7 , - 7?^)], (38) 

where = (F^i-i-Fi 3 )(F^i-i-iAvi)^' stems from the SE diagrams shown in Eigs.|5j)-d, while 77 ^ = (Fi 3 - 7 Avi)^(Fgi-i- 7 Avi)“'(Fg 3 -i- 
/Avi) ' originates from the GSB diagram shown in Eig.|5^. It is interesting to note that the origin of the long-lived oscillations at 
R31, whether predominantly vibrational or vibronic, depends upon the electronic decoherence rates {F^i, Fg 3 ,ri3) and detuning 
Avi = AQ 31 - vi. In Eig. |7^, the ratio \T]el'rig\ between the contributions of the vibronic and vibrational coherences to the 
long-lived beating signal in R31 is displayed as a function of the inter-exciton decoherence rate F 13 , where {F^i, Fg 3 , Avi) are 
taken to be the values estimated from experimental results. Here \rielilg\ > 1 implies that the long-lived beating signal in R31 
is dominated by the vibronic coherence |1) (li| in the electronic excited-state manifold. By fitting the experimentally measured 
beating signals in Nil and R31 to the theoretical model, we found that » 80cm ', which is marked by a vertical dashed 
line in Fig. |7^, where the contribution of the vibronic coherence is ~ 2.5 times greater than the vibrational coherence. These 
results imply that the long-lived beating signal in R31 is dominated by vibronic coherence, originating from electronic excited 
states. It is notable that the vibronic contribution outweighs the vibrational part for a wide range of F^. This is mainly due to 
the fact that the vibronic mixing 62 (*Avi - F 13 ) ' during t 2 depends on the inter-exciton decoherence rate F 13 , while the other 
vibronic mixings ei cx (/Avi h-F^i-F g 3 )“' and 63 cx (-lAvi h-F^i-F^ 3 )“' duringfi and f 3 are independent of F 13 . Considering that 
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FIG. 7. Long-lived beating signals in Nil and R31. a, The ratio Irjelrjg] between the contributions of the vibronic and vibrational coherences 
to the long-lived beating signal in R31. For the experimentally estimated value of Fj 3 , marked by a vertical dashed line, the contribution of the 
vibronic coherence is greater than the vibrational coherence, b, The ratio FjsfFjilTy^ - %!)“' between the amplitudes of the long-lived beating 
signals in N11 (oc F^j“) and R31 (oc F^ 3 F^j*|; 7 j - J 7 j|). In both a and b, we take the values of the parameters estimated from experimental results. 
According to Eq. 1191 . S(Fji -H Fg 3 ) x; 215 cm“* is the theoretical upper bound for F 13 . 


vibronic coherence depends on €2 (see Eq. ([34]i and Figs.| 5]3 and c), while vibrational coherence depends on eie^ (see Eq. ( l36b 
and Eig.|^), the vibronic contribution is increased as ri 3 decreases. We note that these results are in line with the experimental 
observation that the amplitude of the long-lived beating signal in N11 is greater than that of R31 (see Figures 3b and c in the 
main text). In Fig.ITJt, the ratio rg 3 (rgi|? 7 (, - rjgiy^ between the amplitudes of the long-lived beating signals in Nil and R31 is 
displayed as a function of the inter-exciton decoherence rate ri 3 . Here the amplitude of the long-lived beating signal in N11 is 
greater than R31, i.e. rg 3 (rgi| 77 e - > 1, for a range of Eis where the vibronic coherence dominates the long-lived beating 

signal in R31, as shown in Fig.|2t. 


5. Numerical simulation of N11 and R31 

So far the analytic form of the response functions for Nil and R31 were derived with the assumption that the vibronic system 
is well described within the subspace of the vibrational ground and first excited states, which is valid for a small Huang-Rhys 
factor 5 1 . To clarify the validity of this assumption, we performed numerical simulation of the beating signals in N11 and R31 
with higher vibrational excited states, i.e. {|go) , l^i) , ■ ■ ■ , \gn) , |lo), |li) , • • • , |ln) , |3o) , |3i), ■ ■ ■ , |3„)) with n > 1. We found 
that the theoretical beating signals converge for n > 1 and the numerical results are well matched to the analytical results. Here 
the electronic decoherence was modeled by a convex combination of two effective dissipators, i.e. pD\ [p(0] + (1 - 
with 0 < p < 1 , where the dissipators are given by 

A[p(f)] = r,i(2 ID <l|p(f) ID <1| - (ID <1| ,P(0)) + r,3(213> <3|p(f) |3)<3| - {|3> <3| ,p(f))), (39) 

^)2[p(f)] = 2(^11x11+ ^|3><3|)p(f)(^ID<i|+ 7r>l3><3l)-{rgilD<i| + r,3l3><3|,p(f)). (40) 

By substituting electronic coherences |g) {1| and |g) {3| to the dissipators, one can show that both £)i[p(f)] and 2 ) 2 [p(f)] give rise 
to the same set of decoherence rates E^i and r,3 for Ig) {1| and |g) {3|, respectively, implying that the decoherence rates of |g) (1| 
and Ig) (3| are independent of the value of p in the convex combination. For |D {3|, on the other hand, D\ [p(f)] and Diipit)} lead 
to different decoherence rates Fji -i-rg 3 and ( ^jTg\ - Y, respectively. This enables us to vary the inter-exciton decoherence 
rate T 13 within a range of ( ^jTg\ - -y/f^ )^ < ri3 < Fgi - 1 - rg3 by changing the value of p in the convex combination. In addition 
to the electronic decoherence, the relaxation of the vibrational mode was modeled by Eq. (l20l i in the simulations. We found that 
Eq. ( l20b can be approximated by Eq. (1211 1 due to the high vibrational frequency {hv\ » ksT). 


6. Feynman diagrams represented in vibronic eigenbasis 

Here we provide the Feynman diagrams for Nil and R31 represented in the vibronic eigenbasis of the time evolution super¬ 
operator 1/(0, which are equivalent to the Feynman diagrams in the uncoupled state basis shown in Fig.|5] 

For Nil, the vibronic mixing €2 takes place during waiting time t 2 (cf. Fig. [5^), where the vibronic coherences responsible for 
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FIG. 8. Feynman diagrams contributing to the beating signals in Nil and R31 represented in vibronic eigenbasis. a,b, The stimulated 
emission diagrams contributing to the beating signals in N1 1 . c-f, The stimulated emission diagrams contributing to the beating signals in R31. 
g-j, The ground state bleaching diagrams contributing to the beating signals in R31. 


the short-lived and long-lived beating signals in N 11 are given by 

|lo)<3ol = (l+e2V'^"(|lo)<3ol + e 2 |lo)<lil), (41) 

HoXlil = (l+e|r‘^'(|lo><li|-e2|lo)<3ol), (42) 

respectively, where the vibronic eigenstates (3o| oc {3ol + 62 (111 and (li| or (1 1 | - 62 (3o| are normalized by (1 -H not by 

(1 + | 62 p) due to the biorthogonality of the eigenstates of the non-Hermitian operator K in Eq. ( |25]) . When the light-induced 
vibrational excitation of overtones, i.e. {li|, is negligible due to the small Franck-Condon factors, the transition dipole moments 
of the vibronic eigenstates (3o| and (li| are determined by their amplitudes in {3o|, each of which is given by / 23„(1 + ef) and 
respectively. Here denotes the transition dipole moment of {3o|. In the eigenbasis, the Feynman diagrams 
responsible for the short-lived and long-lived beating signals in Nil are described by Figs. [ 8 ^ and b, respectively. Given that 
there are two transitions between (gol and {3o| (and also between (gol and {li|) by the second and third pulses, the square of the 
transition dipole moments of (3o| and (li| is reflected in the response function, each of which is given by/i 3 ^j(l-i-e|)“' jU 3 ^(l-e|) 

and respectively. This is in line with the analytic form of the response function for N11 shown in Eq. OTT l. 

For R31, on the other hand, vibronic mixing takes place during coherence, waiting and rephasing times (fi, f 2 , fa, respectively, 
cf. Figs. |5]r-e). The vibronic mixing ei during coherence time fi leads to the vibronic eigenstates {3 q'| oc (3o| -i- ei (li| and 
(Ij^'l oc (111 - ei {3o|, where vibronic coherences during fi are represented by 

l^o><3'‘*| = (1 +e2)-'/2(|^o)(3o| + e; l^oXhl), 
l^oXll‘^1 = (1 +6?)‘‘^"(l^oXli|-eil^oX3ol). 


(43) 

(44) 
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Here the superindex (1) of (3 q and {1 j*'| reminds us that the vibronic mixing takes place during coherence time ti\ throughout 

~( 2 ) ~( 2 ) 

this work, the vibronic eigenstates {3 q | and {1 j | responsible for the vibronic mixing 62 during waiting time t 2 have, for the sake 
of simplicity, been denoted by (3o| and (1 1 1, respectively. We note that ei in Eq. (IlST i is different from 62 in Eq. (l32T i. as the time 
evolution of the coherences |go) (3ol and |go) (111 during coherence time fi is governed by a different non-Hermitian operator Ki 

Ki = (/n 3 - |3o> <3ol + (/fl, + ivi - r,i - r„) 1 11 > < 11 1 + ivi V^(|3o) < 1 1 1 + 1 1 1 ) <3ol), (45) 

defined by |go) {3ol = Igo) {3o| . In the eigenbasis, the SE diagrams shown in Eigs. |5]r-d can be represented by four 

diagrams shown in Eigs.[ 8 };-f, where the transition dipole moments of (3^'| and are given by /t 3 „(l + and -yU 3 „(l + 

^ 2 )-i/ 2 ^j, respectively. It is notable that the vibronic eigenstates {3 q and during coherence time t\ are different from the 
vibronic eigenstates {3o| and (li| during waiting time t 2 , as the vibronic system is in a superposition between electronic ground 
and excited states (see Eqs. (l43T l and (l44li ') and in the electronic excited-state manifold (see Eqs. (HTt and (|42|i), respectively, 
which leads in general to different values of the vibronic mixings e\ and 62 - The diagrams shown in Eigs.|3;-f describe the fact 
that the vibronic eigenstates (3 q and {lj*^| can be represented by superpositions of (3o| and (Ii|. In Eigs. [ 8 }; and d, for instance, 
the vibronic eigenstate (3 q induced by the first pulse can be represented by a superposition of (3o| and (1 11 

<3''>| = (l+ef)-‘/2((3o| + e,<l,|) (46) 

= (1 + ef)-‘/2(i + e|)-i/2[(i + ejejXSoi + (e, _ 62 )(111)]. (47) 

Here the prefactors of (3o| and (li|, i.e. (1 H-e^)"'^^(l -i-e 2 )”*^^(l +^ 1 ^ 2 ) and (1 -i-e^)“'^^(l H- 62 )~*^^(fi -^ 2 ), enable us to introduce 
two separated diagrams shown in Eigs. | 8 }; and d, where the prefactors are multiplied to the response function, similar to the 
transition dipole moment. Similarly, the other vibronic eigenstate (lj'^| can be represented by a superposition of (3o| and (ii|, 
leading to the prefactors for the diagrams shown in Eigs.| 8 ^ and f. Using the transition dipole moments of (3ol and (li| induced 
by the third pulse, one can show that the response function induced by the SE diagrams is given by Eq. (l34l i: here the lineshape 
functions - Q 3 ) - r^ 3 )“* and (-/(wi - fl 3 -H Avi) - E^O^' along the wi -axis correspond to the diagrams where the vibronic 

system is in |go) <3o’| (c/. Eigs. [ 8 }; and d) and in |go) <1 i'’l (cf. Eigs . | 8 ^ and f), respectively, during coherence time f 1 . 

The vibronic mixing 63 during rephasing time leads to the vibronic eigenstates |3y ) oc |3o) - 1-63 |1 1 ) and |lj ) oc |1 1 )- 63 |3o), 
where vibronic coherences during are represented by 

13®) I = (1 + e3V'^"(|3o> I + 6311 1 ) I), 

Ilf X^ll = (1 + 63 )''^"(| 1 i><^iI - 63 |3o><^ll). 

The time evolution of the coherences |3o) (gi| and |li) (gi| is governed by a non-Hermitian operator 

K 2 = (-*^3 + ivi - r,3 - Tv) |3o> <3ol + - r,i) 1 11 > < 11 1 - ivi V^(|3o) <1 1 1 + 1 1 1 ) <3ol), (50) 

defined by 'Z/(f 3 ) |3o) (gi| = |3o) (gi|. Similar to the SE diagrams, the GSB diagram shown in Eig. can be represented 

by four diagrams shown in Eigs.[ 8 ^-j. Using the transition dipole moments of the vibronic eigenstates, one can show that the 
response function induced by the GSB diagrams is given by Eq. (l3^ . where the lineshape functions (iioji - Qi - Avi) - rg 3 ) ' 
and {i{a >2 - Qi) - E^O^' along the W 3 -axis correspond to the diagrams where the vibronic system is in |3f) (gi| and in |if) (gi|, 
respectively, during tj,. 

These results imply that the Eeynman diagrams for Nil and R31 can be represented in both uncoupled state basis and vibronic 
eigenbasis equivalently, and the analytic form of the response functions in Eqs. (|3TI ). ( l34b and (l36l l is independent of the basis 
chosen to represent the Eeynman diagrams. 


(48) 

(49) 


B. The response function for N22 


Here we provide a vibronic model for bands 2 and 3 of C803, where bands 2 and 3 are coupled to the intramolecular 
vibrational modes with frequency tiV 2 ~ 470 cm“'. 

In Eig. |9^, the absolute square of the Eourier transform of the beating signal in N22 is displayed as a function of a> 2 , which 
is normalized by the amplitude of Nil at Ha )2 ~ 705 cm“'. The amplitude of N22 is maximized around tia )2 ~ 460 cm“* with 
an amplitude in the range of 5 % of the Nil peak. When bands 2 and 3 are coupled to a vibrational mode with frequency V 2 
mediated by an effective Huang-Rhys factor S 2 , the response function for N22 is given by 


Rlgih) ~ ^ilpdn^gl 


^{iAQj2-^23)t2 


+ e 


iV2t2 


/ /V2V^ V 

\/Av 2 - r23 / 


(51) 
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FIG. 9. The relative amplitude of N22 and Nil. a, The absolute square of the Fourier transform of the beating signal in N22 as a function 
of 0 ) 2 , which is normalized to the amplitude of N11 at fioj 2 ~ 705 cm“*. b. Theoretical results of the ratio between N22 and Nil are displayed 
as a function of the Huang-Rhys factor S 2 - The Huang-Rhys factor S i = 0.0006 of the vibrational mode with frequency ftvi s; 668 cm“' is 
marked by a vertical dashed line. 


with Av 2 = AQ 32 - V 2 , l-i 2 p denotes the transition dipole moment of band 2 for light polarized parallel to the longitudinal axis 
of C803 and firg 2 ~ 110 cm ' represents the electronic decoherence rate of band 2, both of which can be estimated using the 
absorption spectrum shown in Fig.| 6 l From the experimentally measured beating signal in N22, we found that /iF 23 ~ 200 cm“^ < 
H^g 2 + Fg 3 ) (not shown). In Fig. |9j), the amplitude of the theoretical N22 is displayed as a function of the Huang-Rhys factor 
S 2 , which is about 5 % of N11 over a range of realistic S 2 values. For a comparison, the Huang-Rhys factor S 1 of the vibrational 
mode with frequency hvi » 668 cm ' is marked by a vertical dashed line. These results imply that the small amplitude of the 
beating signal in N22 is mainly due to the high electronic decoherence rate of band 2. 


C. A correlated fluctuation model for bands 1 and 3 of C803 

Here we provide a correlated fluctuation model for bands 1 and 3 of C803 where coherent interaction between excitons 
and quasi-resonant vibrations is not considered. Within the level of Markovian quantum master equations, we show that the 
experimentally measured long-lived beating signals in N11 and R31 cannot be explained by correlated fluctuations. 

The main idea of the correlated fluctuations is that when bands 1 and 3 are coupled to a common environment, the correlated 
noise enables the inter-exciton coherence |1) {3| to decohere very slowly compared to the coherences |g) (1| and |g) {3| between 
electronic ground state and excitons. This is similar in spirit to the decoherence-free subspaces in quantum information theoryi^. 
Here we consider a Markovian quantum master equation in the form of 

^p(f) = -^[He,p(t)] + X Z - ]^{Al{cS)Ap{cS),p{t )^, (52) 

ijj a,f3 ' ^ 

which is the same to Eq. (3.143) in The Theory of Open Quantum Systems by H.-R Breuer and F. Petruccionei^, which is 
called the Redfield equation with the secular approximation in some literatured^. Here the interaction Hamiltonian is modeled 
by //e_ph = Tia^a ® Ba with Aq, = A^ and Ba = bI, each of which is a Hermitian operator of the system and environmental 
degrees of freedom, respectively. With the exciton states \k) defined by Hg Ik) = SQj. \k), we introduce a projection operator 
n(f2) = \k) {k\ — YjU ^k) \k) {k\ where the Kronecker delta is defined by 6{i, j) = 1 if / = j and 6{i, j) - 0 other¬ 

wise. In other words, n(f2) is a projection operator onto the exciton subspace belonging to the exciton energy Q. In Eq. (l52l l. 
Aa(u) - 2n'-n=w n(Q)AQ.n(Q') = ~ f2)n(Q)AQ,n(Q'). The interaction Hamiltonian N^-ph between excitons and 

background phonons is modeled by Aq. = le^) {ea\ and Ba - 2 ^ hga^(a^^ + af), where ga^ denotes the coupling of the local 
excitation of site a to a background phonon mode f. When ga^ + 0 and g^^ + 0 for different a and /?, spatially separated sites 
a and (3 are coupled to a common phonon mode leading to correlated fluctuations in the energy levels of the different sites a 
and p. The correlated fluctuations are absent when each site is coupled to an independent phonon bath, such that ga^gp^ - 0 for 
all a 4 /? and f. The information of the correlated fluctuations is included in the definition of in Eq. (l52l i 

1 r“ 

j dse‘^UBl(s)Bp(0)}, (53) 

where for fixed a/, japico) form a positive matrisi^. Here yap(<jS) - 0 for all a 4^ if each site is coupled to an independent 
phonon bath and japico) 4 0 for some a 4 y 6 if different sites a and f are coupled to the same phonon modes. 
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Using experimentally measured absorption and 2D spectra of C803, we found that the electronic decoherence rate of the 
coherence |g) {k\ between electronic ground state and band k is given by hFgi 65 cm“' and hFgj 150 cm“' for bands 1 and 
3, respectively. Within the level of the Markovian quantum master equation in Eq. (l52]) . the decoherence rates r,i and rg 3 are 
given by 


r^i = ^ yri^/ + r«i> (54) 

Ml 

= T y 73^/ + 7g3, (55) 

1*3 


where denotes the incoherent population transfer rate from band k to band I, and represents the pure dephasing rate of 
the coherence |g) {k\. The population transfer rates and 73 ^/ can be estimated using the exponential dynamics in 2D spectra. 
To estimate these rates, we performed a global target analysis on all parallel 2D spectra of C803ii. We found that the population 
transfer rates from band 3 to lower energy bands 1 and 2 are approximately given by 73^1 » (300 fs)“* and 73^2 ~ (66 fs)“', 
corresponding to » 18 cm ' and ?J 73-»2 ~ 80 cm respectively, and the other population transfer processes are slow in 

comparison, i.e. yk~>i < (2ps)“'. In this case, the pure dephasing rates of |g) {1| and |g) (3| are given by fiygi » fiFgi » 65 cm“' 
and hygj, ^(rg 3 - 573^1 - 573 -» 2 ) ~ 101 cm“', respectively. 

The electronic decoherence rate T 13 of the inter-exciton coherence 11) {3| between bands 1 and 3 is given by 


ri3 = 



+ 



•/ + 713 , 


(56) 


where 713 is the pure dephasing rate of the inter-exciton coherence in the presence of correlated fluctuations. We found that for 
given 7 gi and 7 ^ 3 , the inter-exciton dephasing rate 713 should be higher than a theoretical lower bound given by 

7i3>(V^-V)^y (57) 

Within the level of the Markovian quantum master equation in Eq. (l5^ . the lower bound is not violated by any spectral densities 
and correlated fluctuations, as will be shown below. Using the estimated values of the pure dephasing rates Hygi ^ 65 cm ' and 
Tiygj, M 101 cm *, we found that the lower bound in Eq. (fSTl) is reduced to S 713 > 4 cm *. Therefore, even in the presence of 
correlated fluctuations, the decoherence rate ri 3 of the inter-exciton coherence |1) {3| in Eq. (l56l) should be higher than a lower 
bound given by 


ri 3 >-(73^1-H73^2)+ () ~(ioofs)'. ( 58 ) 

It is notable that the lowest decoherence rate T^ (lOOfs)^' (c/. tiFi^ » 53 cm“') is too high to explain the long-lived beating 
signals observed in the experiment, as shown in Eig. [TOk . where the simulated results based on the correlated fluctuation model 
are shown as a blue solid line and the experimental results are shown as a light blue line. These results are mainly due to the fast 
population transfer from band 3 to bands 1 and 2 observed in the experiment. 

We note that our results are not sensitive to the estimated values of the population transfer rates. The inter-exciton decoherence 
rate ri 3 is minimized when there is no population transfer between excitons, i.e. y^^i - 0 for all k I, where the coherences 
between electronic ground state and excitons are destroyed only by pure dephasing noise, i.e. fiFgi — tiygi 65 cm ' and 
tiFg-i - Hygi 150 cm“^. Even though this condition is not satisfied for C803, this is the best scenario of the correlated 
fluctuation model where the decoherence rate ri 3 of the inter-exciton coherence is minimized 

ri3=7i3>(V^- Vi^)'^(303fs)-'. (59) 

However, even in this case, the lowest decoherence rate ri 3 (303 fs)~' supported by correlated fluctuations is not low enough 
to explain the experimentally measured long-lived beating signals in N11 and R31, which persist beyond t 2 ~ 800 fs, as shown 
in Eig.fTOb. This is due to the different decoherence rates Fgi < rg 3 of Ig) (1| and Ig) (31 observed in the experiment. This leads to 
a non-zero lower bound on the inter-exciton decoherence rate En, as shown in Eq. (l59l l. In addition, the beating signals in N11 
and R31 consist of a short-lived component with 1/e decay time of ~ 66 fs as well as a long-lived component persisting up to 
f 2 ~ 1 ps. This is contrary to the prediction of the correlated fluctuation model where a single oscillatory component is expected 
with 1/e decay time of For a comparison, the theoretical prediction of the vibronic model is shown in Fig. [TOb . where 
both short-lived and long-lived components are present. We note that in the vibronic model, the decay rate of the long-lived 
component is independent of F^j and Fg 3 , as it is determined by the other degrees of freedom, such as the dissipation rate 7 ,, of 
the vibrations and the degree of vibronic mixing €2 leading to a lifetime borrowing effect dy, as shown in Eq. (l30l) . 
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FIG. 10. Correlated fluctuation model, a. Simulated beating signal for R31 in the presence of correlated fluctuations (without quasi-resonant 
vibrations). The modeled curve is shown as a solid blue line and experimental results are shown as a light blue line. Here we take the lowest 
decoherence rate F^ s; (lOOfs)^' of the inter-exciton coherence allowed by correlated fluctuations, constrained by experimentally determined 
population transfer rates from band 3 to bands 1 and 2. As shown in the inset, correlated fluctuations cannot explain the experimentally 
measured long-lived beating signal in R31. b. Simulated beating signal for R31 in the absence of the exciton relaxation, shown as a solid 
blue line. The correlated fluctuation model predicts the lowest decoherence rate F^ s; (303fs)“* of the inter-exciton coherence when there 
is no exciton relaxation, even though this condition is not satisfied for C803. Nevertheless, the correlated fluctuation model cannot explain 
the experimentally measured long-lived beating signal in R31, which persist beyond t 2 ~ 800 fs, as shown in an inset, c. Simulated beating 
signal for R31 in the presence of a vibrational mode with frequency hvi ~ 668 cm“', which is quasi-resonant with the exciton energy splitting 
AHsi between bands 1 and 3. Here the vibronic and vibrational coherences induce a long-lived beating signal in good agreement with the 
experimental results. The root-mean-square deviation (RMSD) of the experimental results and theoretical prediction in a, b and c is 0.74, 0.86 
and 0.59, respectively. We note that the correlated fluctuation model can also not explain the long-lasting beating signal in N11 (not shown). 


We now derive Eqs. (I54b-(l57b using the Markovian quantum master equation in Eq. ( l52b . 

(Dephasing noise) We start with the case that w = 0, leading to the pure dephasing noise. Eor the sake of simplicity, we assume 
that there is no degeneracy in the exciton energies Q.ic, such that + Q; for all k ^4 /. In this case, Aq,( 0) = Xit 1^) (^1 1^) (^1 = 
Fik K^ka)P \k) {k\. By substituting |g) <1|, |g) (3| and |1) (3| to the dissipator of the quantum master equation for to - 0 

^ r.M0) (A/5(0)p(f)Al(0) - i{A^(0)A^(0),p(f))), (60) 

one obtains the following pure dephasing rates of the coherences |g) {1|, |g) {3| and |1) {3| 

^ X l<3k«)lV«M0) K3k/^>l' - |V3|", 

afi 

ri3 = ^ " meaf)yai3{0)[\{l\ept - K3k/5)l") = |vi - V3l' , 

afi 

where v 4 represents a vector dehned by v 4 = here is a real vector with elements Kkka)|^ > 0 and y is a positive 

matrij&i^ with elements -/a/iiO), leading to to a positive matrix dehned by y = yh 2 .^i /2 Pqj- gjyeji puj-g dephasing rates 
2 2 

7gi — kil and yg 3 = |y 3 | , the inter-exciton pure dephasing rate y^ is constrained by 

yi3 = |vi - V3p > (|vil - |V3l)^ = ( V^- f ’ (64) 

due to the triangle inequality, |vi - V 3 I -H |v 3 | > |vi| and |vi - V 3 I + |vi| > |v 3 |. Here the equality holds if and only if vj is parallel 
to V 3 , which depends on y*^^ (spectral densities and correlated fluctuations) as well as wi and VV 3 (the spatial overlap between 
excitonic wavefunctions). Note that the lower bound in Eq. dMl i has been derived based on the positivity of y'^^, which is 
satisfied for any spectral densities and correlated fluctuations. These results imply that the inter-exciton dephasing rate yi 3 can 
be reduced by the correlated fluctuations as well as the spatial overlap between excitonic wavefunctions. 

(Exciton relaxation) We now consider the case that w it 0, leading to the incoherent population transfer between excitons. 
With AQ.ici - Qk - Q.I denoting the exciton energy splitting between bands k and /, Aq(w) = 2i/t,/<5(w, AQ^./) |() {/| Aq |k) {k| = 
Yjkj <5(w, AQ.ki) (/ka) {ea\k) \l) {k\ and the dissipator of the quantum master equation for w it 0 is given by 

^p(f) = zi; japlAElki) i^p{AQ.ki)p{t)A\{AQ.ki) - -{A^(AQH)A^(AQH),p(f))j , 


( 61 ) 

(62) 

(63) 


(65) 
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where the population transfer rate from band k to band / is given by 

yk^i = ^ {k\ea) {ea\l) Jaisi^^ki) {l\e/j){e/j\k) > 0, (66) 

a,/3 

which is positive, as 7 q,^(AQ^/) form a positive matri?&i2 for given AQ/./. By substituting |g) (1|, \g) (3| and |1) (3| to the dissipators 
in Eqs. ( l60b and < l65] t, one can show that the electronic decoherence rates Fgi, rg 3 and ri 3 satisfy Eqs. (I54lt-<l57]t. These results 
are valid for any spectral densities and correlated fluctuations within the level of the Markovian quantum master equation in 
Eq. (l52l t, which includes the theoretical models considered in previous studies^. 

In summary, the correlated fluctuation model cannot explain the long-lived beating signals in N11 and R31 within the level 
of Markovian quantum master equations, as the decoherence rate ri 3 of the inter-exciton coherence 11) (3| is constrained by the 
experimentally observed asymmetric decoherence rates Egi and rg 3 , i.e. Tg^ ^ ^Egi, and the fast population transfer from band 
3 to bands 1 and 2. We note that the asymmetric decoherence rates Tgi and rg 3 are related to the fact that i) the lineshape of R31 
is elongated along a»i-axis (cf. Eigures 2b and d in the main text), ii) the amplitudes of the short-lived beating signals in Nil 
and R31 are different in magnitude (cf. Eigures 3b and c in the main text), and iii) the absorptive linewidths of bands 1 and 3 are 
different (cf Eig.[ 6 ]in the SI). 
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